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FOREWORD 

~ —  ^The  present  report  is  part  of  a  two  volume  set  which  describes 
a  nonlinear  solid  rocket  motor  instability  analysis  and  computer  pro¬ 
gram.  Volume  I  contains  the  analytical  basis  for  the  computer  program 
and  a  discussion  of  the  results  obtained  to  date:  Volume  II  of  the  set 
describes  the  computer  program  and  serves  as  a  user's  manual. 

This  investigation  is  entitled  NUMERICAL  ANALYSIS  OF 
NONLINEAR  LONGITUDINAL  COMBUSTION  INSTABILITY  IN  METALIZED 
PROPELLANT  SOLID  ROCKET  MOTORS.  The  two  volumes  are  additionally 
subtitled  as  follows: 

Volume  I  -  Analysis  and  Results 

Volume  II  -  Computer  Program  User's  Manual 

This  investigation  was  sponsored  by  the 

Air  Force  Rocket  Propulsion  Laboratory 

Director  of  Laboratories 

Edwards,  California  93523 

Air  Force  Systems  Command ,  United  States  Air  Force 
under  contract  number  F04611-71-C-0060  with  Robert  J.  Schoner  as 
technical  monitor.  Jay  N.  Levine  of  Ultrasystems  (formerly  Dynamic 
Science)  was  program  manager. 

This  technical  report  has  been  revievv’ed  and  is  approved. 


Paul  J.  Daily,  Lt.  Col.  USAF 
Chief,  Technology  Division 
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ABSTRACT 


Ihe  primary  objective  of  the  current  effort  was  the  development  and 
SO'  Ion  of  a  nonlinear  analytical  longitudinal  instability  model,  which 
would  allow  all  of  the  various  governing  phenomena  to  be  accounted  for  in 
a  coupled  manner.  The  two  primary  elements  of  the  current  instability 
analy&iS  are  a  method  of  characcaristics  solution  of  the  two  phase  flow 
in  the  combustion  chamber  of  the  motor,  and  a  coupled  calculation  of  a 
transient  burning  rate.  The  transient  burning  rate  analysis  presented, 
herein,  is  a  unl-^ue  and  interesting  development.  It  is  based  on  an  extension 
rf  che  most  popular,  linear,  harmonic  combustion  response  model.  The 
o  r.rent  method  e -lows  the  calculation  of  propellant  burning  response  to  a 

i 

pressure  disturbance  of  arbitrary  waveform,  for  all  time,  including  the 
period  immediately  fc  llowing  the  Initiation  of  the  disturbance.  The  analysis 
also  Includes  a  model  for  velocity  coupled  response.  Therefore,  for  the  first 
time,  the  nonlinear  effects  of  velocity  coupling  on  the  growth  of  pressure 
waves  in  a  combustion  chamber  can  be  computed . 

The  instability  solution,  itself,  begins  with  the  calculation  of  the 
steady  state  two-phase  flow  in  the  motor.  The  flow  in  the  combustion 
chamber  is  calculated  by  numerically  integrating  the  equations  of  motion. 

The  nozzle  flow  is  found  using  the  constant  fractional  lag  approximation.  The 
steady  state  conditions  are  then  perturbed  and  the  subsequent  wave  motion  in 
the  motor  is  calculated  numerically,  using  the  method  of  characteristics. 

The  nature  of  the  engine  response  is  dependent  up  n  the  interaction  the 
various  gain  and  loss  mechanisms  in  the  engine,  which  are,  in  turn,  a 
function  of  the  propellant  burning  response,  the  size  and  amount  of  parti¬ 
culate  matter  present,  the  magnitude  and  shape  of  the  initial  disturbance 
and  the  geometrical  configuration  of  the  motor. 

The  instability  model  is  currently  subject  to  the  following  limitations. 
Only  motois  with  cylindrlcally  perforated  grain  were  considered.  The  gas- 
dynamic  flow  was  assumed  to  be  one-dimensional  and  the  particles  in  the 
gas  stream  were  taken  to  be  of  uniform  size  and  inert.  The  nozzle  flow  is 
assumed  to  be  quasi-steady. 
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A  series  of  instability  solutions  have  been  calculated,  wherein  some 
of  the  main  parameters  such  as  particle  size,  burning  rate  constants,  and 
initial  disturbance  waveform  and  magnitude  have  been  varied,  in  an  attempt 
to  qualitatively  assess  the  behavior  and  validity  of  the  present  model. 

From  all  appearances ,  the  behavior  of  the  model  is  quite  realistic  and 
limited  comparisons  with  data  have  been  quite  encouraging. 
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NOMENCLATURE 

-  burning  rate  parameter,  Eq.  (4-14) 

-  admittance  function 

-  nozzle  admittance  function 

-  gas  only,  sound  speed 

-  sound  speed  based  on  Pp  and  Tp 

-  burning  rate  parameter,  Eq.  (4-32) 
also,  fractional  lag  parameter,  Eq.  (5-6) 

-  burning  rate  parameter  for  velocity  coupling 

-  ratio  of  solid  to  gas  specific  heats,  C  /C 

s  p 

-  constant  in  steady  state  burning  rate,  Eq.  (3-5) 

-  particle  drag  coefficient 

-  erosive  burning  constant,  Eq.  (3-5) 

-  see  cij 

-  specific  heat  of  gas  at  constant  pressure 

-  specific  heat  of  solid  particles 

-  specific  heat  of  gas  at  constant  volume 

-  defined  in  Eq.  (7-34) 

-  port  diameter 

-  normalized  surface  activation  energy,  E  /R  T 
also,  integral  defined  by  Eq.  (8-16) 

-  activation  energy  of  surface  reaction 

-  internal  energy 

-  particle-gas  interaction  force  per  unit  volume,  Eq.  (3-8) 

-  frequency 

also,  as  defined  by  Eq.  (8-14) 

-  defined  by  Eq.  (7-8) 

-  defined  by  Eq.  (6-17) 

-  defined  by  Eq.  (4-20) 

-  enthalpy 

-  defined  by  Eq.  (8-12) 

-  fractional  lag  constant,  Eq.  (5-3) 

-  chamber  fractional  lag  constant,  Eq.(7-44) 

-  thermal  conductivity,  also  complex  wave  number 
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thermal  conductiviti'  of  the  solid  particles 

length  of  the  grain,  also  fractional  lag  constant,  Kq.  (5-4) 

chamber  fractional  lag  constant,  Eq.  (7-44) 

perimeter  of  the  grain 

Mach  number,  also  number  of  points  on  initial  line 
Mach  number  at  burning  surface 
particle  mass,  also  surface  mass  flux 
Nusselt  number 

pressure  exponent  in  steady  state  burning  rate 
constant  in  velocity  coupled  analysis,  Eq.  (4-75) 
exponent  of  pressure  dependence  of  surface  reaction  rate 
pressure 

reference  pressure  in  steady  state  burning  rate 
chamber  pressure 
Prandtl  number 

pressure,  also  used  for  Laplace  transform  variable 
defined  by  Equation  (8-10) 
heat  release  per  unit  mass 

particle-gas  heat  transfer  rate  per  unit  volume,  Eq.  (3-14) 
heat  of  reaction  for  processes  at  burning  surface 
see  ; 

gas  constant,  also  normalized  throat  radius  of  curvature 
universal  gas  constant 
response  function ,  Eq.  (4-35) 

Reynolds  number  based  on  particle  diameter  and  particle-gas 
relative  velocity 

right  hand  side  of  a  characteristics  compatibility  relation 
linear  burning  rate 
area  of  burning  surface 

dimensionless  Laplace  transform  variable,  =  iH  j'/r‘ 
temperature 

adiabatic  flame  temperature 
time 

defined  by  Eq.  (7-33) 
axial  velocity 
threshold  velocity 
defined  by  Equation  (4-28) 


w 

X 

a 


8l 

^2 


Y 

6 

6’ 


®1'  ®2 
)i 

s 

A 

X 

a 

"s 

0 


0 

s 

a 


Cl,  C2 


Tt 

$ 

cp 

0 

uu 


-  reaction  rate  divided  by  gas  density 

-  axial  distance 

-  growth  constant 

-  particle  damping  constant 
“  defined  in  Eq.  (8-35) 

-  particle  to  gas  weight  flow  ratio 

-  ratio  of  particle  to  gas  mass  burning  rates,  Wp/w 

-  ratio  ol  ^as  specific  heats, 

-  a  small  increment  in  time 

-  equal  to  t  6 

-  convergence  criteria  for  characteristics  calculations,  also  used 
in  velocity  coupling  analysis  (Eq.  4-72) 

-  thermal  diffusivity  of  the  propellant 

-  defined  by  Eq.  (4-27) 

-  complex  function  of  frequency,  Eq.  (4-8) 

-  viscosity 

-  equals  fxA 

s 

-  density 

-  density  based  on  Pp  and  Tp 

-  density  of  the  metal  oxide  particles 

-  density  of  the  solid  propellant 

-  particle  radius 

-  defined  by  Eq.  (4-62) 

-  nondlmensional  time,  r®t/4H  also  used  in  Section  2  to 
denote  period  of  oscillation 

-  characteristic  relaxation  time  for  particle  velocity,  Eq.  (3-1) 

-  characteristic  relaxation  time  for  particle  temperature,  Eq.  (3-1) 

-  defined  by  Eq.  (5-10) 

-  phase  angle 

-  nondlmensional  frequency,  Eq.  (4-9) 

-  mass  burning  rate,  per  unit  length,  per  unit  cross-sectional 
area,  Eq.  (3-4);  also  occasionally  used  for  angular  frequency 
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Subscripts 
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end  of  the  propellant  grain 

flame 

gas 

for  the  £th  mode  of  oscillation 
particle 

initial  or  stagnation  value 
at  the  nozzle  throat 

at  the  burning  surface  of  the  propellant 


Superscripts 


( )*  -  in  Sections  3,  5  and  6  only,  denotes  a  dimensional  '*j:iable 

(  )'  -  denotes  fluctuation 


() 


()“ 

(,(r) 


in  Section  4  denotes  steady  state  variable,  in  Section  5 
denotes  an  "equivalent"  gas  value,  in  Section  7  denotes 
an  average  quantity 

pertaining  to  a  right  running  characteristic 
pertaining  to  a  left  nmning  characteristic 
real  part  of 
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1. 


INTRODUCTION 


Solid  propellant  rocket  motors  are  often  subject  to  chamber  pressure 
oscillation  during  the  course  of  their  firing.  The  existence  of  these 
fluctuations  is  indicative  of  the  existence  of  a  coupling  mechanism,  or 
mechanisms,  between  the  combustion  at  the  surface  of  the  propellant  and 
the  gasd/namic  flow  in  the  chamber.  When  the  rate  at  which  energy  is 
supplied  to  the  flow,  by  this  coupling,  exceeds  the  rate  at  which  energy 
is  lost  through  the  various  dissipative  mechanisms  that  exist,  the  chamber 
pressure  oscillations  are  amplified.  They  may  grow  to  such  amplitudes 
that  the  oscillating  acceleration  of  the  vehicle  may  produce  failure  of  the 
equipment  or  even  failure  of  the  motor. 

Of  the  various  unstable  motions  observed  in  solid  propellant  rocket 
chambers,  longitudinal  combustion  instability  is  currently  the  most 
troublesome.  The  axial,  or  longitudinal,  mode  of  combustion  Instability 
•  occurs  in  the  frequency  range  Intermediate  between  the  very  low  frequencies 
of  bulk,  or  L*,  instabilities,  and  the  high  frequencies  associated  with  the 
tangential,  or  transverse,  instability  modes.  Unlike  transverse  mode 
instabilities,  which  are  very  often  eliminated  by  the  addition  of  powdered 
aluminum  to  propellant  formulations;  the  existence  of  a  remedy  of 
comparable  simplicity  for  longitudinal  instabilities  has  not  yet  been 
demonstrated. 

Efforts  to  deal  with  the  problem  of  longitudinal  Instability,  most 
logically,  have  begun  with  experimental,  and  analytical,  investigations 
designed  to  shed  light  on  the  governing  physical  phenomena.  Both 
laboratory  scale  experiments,  and  full  scale  firings,  have  identified  the 
existence  of  numerous,  complex,  mostly  nonlinear,  processes  which  actively 
play  a  role  in  determining  the  stability  characteristics  of  a  rocket  motor. 

As  a  result  of  this  complexity,  investigators  attempting  to  predict  these 
disturbances  have,  in  the  past,  been  almost  universally  forced  to  linearize 
the  equations  governing  the  various  phenomena.  Attempts  to  solve  more 
realistic  nonlinear  models  have  been  made,  but,  only  by  focusing  on  one 
of  the  many  processes  at  a  time.  In  view  of  this,  the  primary  objective  of 
this  program  was  the  development  and  solution  of  a  nonlinear  analytical 
model  for  describing  longitudinal  mode  instability,  which  would  allow  all 
the  various  processes  to  be  accounted  for  in  a  coupled  manner.  A  second 
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analyses  of  the  burning  rate  response  to  harmonic  disturbances. 

The  formulation  of  the  boundary  and  initial  conditions  for  the  problem  is 
discussed  in  Sections  5  and  6.  The  numerical  methods  employed  in  solving 
the  governing  equations  are  discussed  in  Section  7;  while  in  Section  8, 
a  linear  stability  analysis  is  developed,  for  the  purpose  of  comparison 
with  the  nonlinear  numerical  results. 

Time  has  limited  the  number  of  instability  solutions  which  could  be 
obtained  during  the  course  of  the  present  investigation.  However,  as 
shown  in  Section  9,  the  calculated  results  are  quite  interesting,  and  exhibit 
all  the  proper  qualitative  trends.  The  nonlinear  analysis  described  herein 
appears  to  have  the  potential  to  lead  the  way  to  greater  understanding  of 
longitudinal  combustion  instability,  and  to  the  realization  of  more  accurate 
quantitative  predictive  capability. 

A  further  discussion  of  the  conclusions  that  can  be  drawn  from  the 
present  study  may  be  found  in  Section  10.  At  the  end  of  the  report,  in 
Appendix  A,  the  behavior  of  the  transient  burning  rate  model  is  examined 
for  some  special  cases. 


2. 


BACKGROUND  AND  SURVEY  OF  EXISTING  EXPERIMENTAL  AND 
ANALYTICAL  RESULTS 


The  problem  of  acoustic  waves  in  a  closed  tube  is  of  course  a  very 
old  one.  Only  those  recent  works  concerned  with  finite  amplitudes  are 
cited  here.  In  simplest  form,  the  experimental  apparatus  consists  of  a 
rigid  tube  closed  at  one  end  and  fitted  with  a  piston  at  the  other.  All 
reported  experimental  results ,  and  all  but  one  of  the  analyses  give  information 
only  about  the  waves  after  limiting  amplitude  has  been  reached. 


A  periodic  motion  can  be  excited  in  the  tube  at  any  frequency,  but  the 
most  interesting  behavior  occurs  when  the  piston  oscillates  with  a  frequency 
at  or  near  one  of  the  resonant  frequencies  of  the  tube.  It  has  long  been 
known  (see  Ref.  1  and  other  v,;orks  cited  there)  that  as  the  amplitude  is 
increased,  the  wave  motion  in  the  tube  changes  from  a  simple  sinusoidal 
motion  in  both  space  and  time,  into  one  having  both  a  continuous  and  a 
discontintious  part.  At  sufficiently  high  amplitude,  one  or  more  weak  shock 
waves  propagate  to  and  fro.  The  appearance  of  weak  "discontinuities"  is  a 
consequence  of  the  nonlinear  convective  effects  and  the  dependence  of  the 
speed  of  sound  on  temperature. 


Analyses  of  the  steady  wave  motion  as  a  weak  shock  embedded  in 
a  continuous  periodic  motion  have  been  given  by  Saenger  and  Hudson  (Ref.  1) 
and  Betchov  (Ref.  2).  Subsequently,  Chester  (Ref.  3)  discussed  in  more 
detail  the  relative  importance  of  the  mechanisms  involved  and  produced  a 
single  solution  valid  over  a  range  of  frequencies  near  resonance.  He  showed, 
in  agreement  with  Betchov,  that  the  effects  of  heat  conduction  and  viscous 
forces  provide  small  corrections  to  the  Inviscid  solution  but  are  not  required 
to  limit  the  amplitude  of  the  motion. 


In  Ref.  4,  Coppens  and  Sanders  analyzed  the  motion  in  terms  of  the 
generation  of  harmonics.  Their  results  are  not  carried  far  enough,  explicitly, 
to  apply  to  a  discontinuity.  They  also  quote  measurements  which,  as  do 
other  experimental  observations,  show  an  important  feature,  namely  the 
appearance  of  harmonic  distortion  of  the  waverform  at  quite  modest  ampli¬ 
tudes.  For  example,  when  the  amplitude  of  the  fundamental  mode  is  approxi¬ 
mately  one  percent  of  the  ambient  pressure,  the  amplitude  of  the  second 
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harmoric  is  already  fifteen  percent  of  the  amplitude  of  the  first  harmonic. 

A  discontinuity  is  apparent  when  the  amplitude  is  one-tenth  of  the  ambient 
pressure.  This  is  a  general  result  for  acoustic  waves  in  resonant  tubes — 
nonlinear  behavior  is  obvious  at  these  low  amplitudes. 

On  the  other  hand,  in  solid  propellant  motors,  quite  clean  sinusoidal 
motions,  with  very  little  harmonic  content,  are  often  observed  to  amplitudes 
as  high  as  20%  or  more  of  the  average  pressure.  The  explanation  for  this 
has  not  been  definitely  established,  although  it  mus*  evidently  rest  on  the 
behavior  of  the  combustion  processes. 

The  stability  of  wave  motions  in  a  rocket  motor  is  a  primary  problem. 

It  is  then  necessary  to  examine  the  transient  growth  of  waves  from  some 
specified  initial  disturbance.  The  corresponding  problem  for  waves  in  a 
resonant  tube  has  been  treated  in  only  one  work.  Ref.  5.  Perhaps  the  main 
reason  for  this  is  that  the  situation  arises  only  for  a  self-excited  oscillating 
system,  and  hence  is  not  of  interest  when  waves  are  driven  by  external  means. 
The  computation  of  Ref.  5  does  not  include  the  influence  of  combustion  or 
mean  flow,  and  hence  shows  only  the  transient  development  of  the  weak 
discontinuity.  Numerical  results  were  reported  for  a  finite-difference 
calculation  and  for  the  method  of  characteristics. 

As  an  aid  to  interpreting  later  results  for  the  behavior  of  nonlinear  waves 
in  a  motor,  it  is  useful  to  sketch  an  idealization  of  the  motions  in  a  resonant 
tube.  Measurements  of  pressure  are  usually  made  at  the  end  of  a  chamber, 
but  it  is  also  interesting  to  see  the  distribution  of  pressure  along  the  chamber 
as  it  varies  in  time.  Three  cases  will  be  examine:  a  purely  sinusoidal 
oscillation  at  the  fundamental  frr  uency,  a  week  shock,  and  the  sum  of 
these  two.  These  are  shown  in  Figures  2-1  to  2-3,  with  the  distribution 
of  pressure  shown  for  every  eighth  cycle.  The  position  of  the  shock 
wave  relative  to  the  sinusoidal  oscillation  has  been  chosen  such  that  it 
arrives  at  one  end  just  as  the  sinusoidal  pressure  at  that  end  is  a 
maximum.  This  will  not  in  general  be  true  in  actual  cases,  owing  to 
the  influence  of  energy  losses,  average  flow  and,  of  course,  strongly 
nonlinear  effects.  Also,  for  the  example  shown  here,  the  shock  is  assumed 
to  be  sufficiently  weak  that  it  propagates  always  at  constant  speed  equal  to 
the  sound  speed  for  the  standing  wave;  the  period  of  both  motions  is  T. 
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Figure  2-3.  Distribution  of  Pressure  Along  the  Chamber  and  Pressure 
at  X  =  L  for  a  standing  wave  with  a  weak  travelling 
shock  wave 


Note  that  if  the  standing  wave  is  weak,  so  the  motion  is  dominated  by  a 
travelling  shock  wave,  the  pressure  measured  at  one  end  appears  roughly 
as  a  sawtooth  or  triangular  wave.  In  time,  the  pressure  jumps  abruptly,  and 
then  decays  until  it  jumps  again  upon  the  next  arrival  of  the  shock.  Many  of 
the  pressure  traces  reported  in  the  work  at  GARDE  and  SRI,  discussed  below, 
exhibit,  qualitatively,  this  kind  of  behavior.  Even  for  the  case  of  nonlinear 
motions  in  a  simple  resonant  tube,  the  pressures  measured  differ  from  the 
form  shown  here,  owing  to  contributions  from  many  harmonics. 
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It  is  clear  that  previous  work  on  the  classical  problem  of  waves  in  a 
resonant  tube  provides  background,  but  no  useful  results  for  the  problem  of 
transient  waves  in  a  rocket  motor. 


2 . 2  Experimental  Results  for  Longitudinal  Waves  in  Solid  Propellant 

Rocket  Motors 

Longitudinal  instabilities  arise  spontaneously  in  full-scale  motors — 
a  number  of  examples  exist,  but  few  have  been  reported  in  the  open  literature. 
The  most  recent  reports  of  such  observations  appear  in  Refs.  6-8.  Most  of  the 
data  which  are  available  for  motors  have  been  taken  in  laboratory  devices, 
and  almost  all  involve  pulsing. 

It  is  possible  to  produce  sharp  fronts  propagating  in  a  tube  by  introducing 
an  impulse  of  mass,  momentum,  or  energy.  This  may  be  done,  for  example, 
with  an  electric  spark  or  a  small  explosive  charge.  If  the  tube  contains  only 
an  inert  gas,  the  subsequent  motions  are  not  particularly  interesting — they 
die.  For  a  rocket  motor,  however,  introducing  a  pulse  during  a  firing  is  a 
useful — indeed,  the  only — means  for  determining  completely  the  stability 
characteristics.  The  combustion  processes  constitute  a  source  of  energy; 
the  amplitude  of  the  pulse  may  grow,  remain  unchanged,  or  decay.  Whatever 
happens  must,  of  course,  reflect  the  nonlinear  conflict  between  the  gain  and 
loss  of  ene  gy  for  the  pulse. 


Qualification  of  liquid  rocket  motors  by  pulsing  has  been  a  standard 
procedure  for  some  time.  Although  the  technique  has,  for  practical  purposes 
essentially  not  been  used  for  solid  propellant  motors  (an  exception  is  program 
for  developing  motors  used  in  the  Canadian  "Black  Brant"  vehicle,  cited  in 
Ref.  17),  there  exists  a  substantial  amount  of  laboratory  data  for  instability 
studies  using  T-burners  as  well  as  burners  with  cylindrical  configurations. 
Only  results  for  end- vented  burners  will  be  discussed  here. 

It  should  be  noted  that  in  laboratory  tests,  the  strength  of  the  pulse 
can  be  controlled,  although  precision  is  a  problem.  Thus,  it  is  possible 
to  produce  perturbations  much  larger  than  one  would  normally  encounter  in 
an  actual  motor,  as,  for  example,  might  occur  due  to  partial  blockage  of  the 
nozzle  by  a  small  piece  of  material.  Laboratory  tests  therefore  offer  more 
flexibility,  and  by  proper  interpretation  should  be  applicable  to  real  motors. 
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A  significant  amount  of  experimental  data  has  been  taken  at  the 
Canadian  Armament  and  Research  Institute  (GARDE)  (see  Refs.  9-18)  and 
at  the  Stanford  Research  Institute  (Refs.  19-22).  More  recently,  work  on 
small  pulsed  motors  has  been  done  at  the  Aerojet-General  Solid  Propulsion 
Company  (Ref.  23).  Limited  results  on  spontaneous  longitudinal  oscillations 
have  been  obtained  at  the  Naval  Weapons  Center  (Ref.  24),  for  a  center- 
■  vented  configuration .  A  few  observations  were  reported  earlier  by  Hercules , 
Inc.  (Ref.  25)  and  by  the  Ballistics  Research  Laboratory  (Ref.  26). 

2.2.1  Data  Taken  at  CARDE  and  SRI 

In  both  the  CARDE  and  SRI  work,  pulses  were  introduced  at  the  head  end 
of  a  motor  by  introducing  small  explosive  charges.  Apart  from  a  few  important 
visual  observations  (Refs.  13,  15  and  16)  the  data  consists  of  pressure 
measurements  taken  at  the  head  and  aft  ends  of  the  chamber.  Rectangular, 
slotted,  star,  and  other  cross-sections  have  been  used,  but  most  of  the 
tests  have  been  made  with  circular  ports.  Experimental  variables  which 
have  been  studied  include  chamber  diameter  and  length  (hence  frequency), 
throat/port  area  ratio  (hence  chomber  pressure  and  Mach  number),  initial 
temperature,  and  composition  of  the  propellants.  The  last  involved  ballistic 
additives  and  aluminum  content  as  well  as  changes  of  binder  and  size  and  type 
of  oxidizer. 

Pressure  traces  suggest  ihe  presence  of  one  or  more  discrete  finite 
pressure  waves  propagating  up  and  down  the  chamber.  (See Figures  1-3  and 
following  remarks) .  Subsequently,  optical  observations  (Refs.  13,  15,  and 
16)  established  that  indeed  the  instability  involved  a  weak  shock  wave, 
usually  travelling  with  Mach  number  less  than  1.2.  It  is  evident  from  the 
experimental  results  that  the  disturbance  introduced  by  the  explosive  charge 
very  rapidly  settles  to  the  nearly  planar  wave  travelling  to  and  fro  in  the 
chamber.  The  strength  is  a  function  of  position  in  the  chamber,  mainly 
because  of  the  losses  incurred  upon  reflection  at  the  ends,  particularly 
at  the  nozzle.  The  influence  of  the  mean  flow  was  more  pronounced  during 
travel  from  the  nozzle  to  the  head  end;  the  wave  front  was  then  concave  in  the 
direction  of  propagation. 


The  wave  motion  therefore  appears  to  be  analogous  to  the  case  of 
the  classical  resonant  tube  driven  at  high  amplitude,  an  idealized  example 
of  which  is  shown  above  in  Figures  2-1  to  2-3.  However,  in  a  motor,  there 
are  obviously  strong  perturbations  due  to  the  combustion  and  related 
processes.  That  the  instability  is  often  predominantly  a  discrete  travelling 
wave  (in  some  cases  two  waves  are  present)  is  probably  related  to  the 
pulsing,  but  insufficient  information  is  given  in  the  accounts  to  draw  any 
conclusions.  In  the  GARDE  work,  the  explosive  charges  used  were  0.3  gm 
to  12  gm  depending  on  the  size  of  the  motor.  No  further  details  are  given. 

Moreover,  only  portions  of  a  few  records  are  reproduced  and  all  show 
the  discrete  waves  very  shortly  after  the  pulse.  It  is  not  possible  to  deduce 
any  information  about  growth  or  decay  rates  of  the  waves.  Indeed,  in  the 
latest  work  reported  (Ref.  18),  only  time  averaged  pressures  were  recorded. 

The  experimental  results  therefore  relate  to  the  qualitative  question 
of  whether  a  pulse  is  stable  or  unstable,  although  in  a  few  cases  the  limiting 
amplitudes  are  given.  The  severity  of  the  instability  is  interpreted  usually 
in  terrfis  of  the  shift  of  mean  pressure. 

In  the  motors  used,  progressive  burning  causes  the  mean  pressure  to 
rise  dOring  a  firing;  testing  the  stability  characteristics  consists  of  intro¬ 
ducing  a  sequence  of  pulses  which  therefore  occur  at  successively  higher 
pressures.  Typically,  pulses  early  in  a  firing,  at  lower  mean  pressures, 
may  be  stable,  while  those  at  the  higher  pressures  late  in  the  firing  may 
produce  the  travelling  wave  instability.  (It  is  of  course  possible  that  in 
some  cases  all  pulses  may  be  stable  or  all  unstable).  Since  there  is 
generally  an  increase  of  the  mean  pressure  when  the  pulses  are  unstable, 
there  are  some  delicate  questions  involved  in  determining  a  stability 
boundary  with  any  precision.  How  this  has  been  done  in  the  work  cited 
(especially  Refs.  13  and  18)  will  not  be  discussed  here;  it  is  thoroughly 
discussed  in  those  references) . 
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The  behavior  of  the  mean  pressure  during  a  firing  in  which  an  unstably 
pulse  is  initiated  is  sketched  in  Figure  2-4 
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Figure  2-4.  Mean  Pressure  for  a  Typical  Firing  from 
the  GARDE  and  SRI  Work 


If  no  pulse  is  introduced,  the  firing  is  normally  stable;  thus  the  data  truly 
relates  to  a  problem  of  nonlinear  instability.  Much  of  the  interpretation  of 
data  is  based  on  plots  of  mean  chamber  pressure,  measured  at  a  fixed  fraction 
of  web  burned,  versus  K  ,  the  ratio  of  burning  surface  to  throat  area.  An 
example  is  shown  in  Figure  2-5;  such  data  can  be  obtained,  for  example, 
from  a  series  of  firings  identical  except  for  the  initial  value  of  throat  to 
port  area  ratio. 


Locus  of  unstable 
operation 


Mean 

Chamber 

Pressure 


Per 


K  _  Burning  Surface  Area 
n  Nozzle  Throat  Area 


Figure  2-5.  Sketch  of  Stable  and  Unstable  Operating 
Pressures  for  a  Fixed  Fraction  of  Web 
Consumed 


The  intersection  of  the  two  lines  in  Figure  2-5  gives  "critical"  or  "threshold" 
values  of  pressure  and  K^. 

Early  work  seemed  to  suggest  that  the  dynamical  characteristics  of  many 
propellants  could  be  correlated  in  terms  of  the  threshold  values  of  pressure 
and  (Refs.  10,12  and  21),  However,  subsequently  (Ref.  18)  it  was 
established  that  not  only  was  this  correlation  severely  limited  by  the  motor 
configuration,  but  also  that  by  no  means  all  propellants  could  be  characterized 
by  the  same  stability  boundary.  Even  worse,  in  many  cases  the  stable  ana 
unstable  lines  shown  in  Figure  5  had  slopes  so  nearly  alike  that  a  well-defined 
threshold  point  could  not  be  established.  As  a  result  of  the  last  conclusion, 
Roberts  and  Brownlee  (Ref.  18)  were  forced  to  interpret  their  data  in  a 
different  way.  They  chose  to  use  the  difference  of  the  mean  pressures  for 
stable  and  unstable  operation,  taken  at  some  arbitrary  value  of  the  stable 
value  (1200  psia). 


It  is  therefore  not  possible  at  the  present  time  to  discuss  this  data 
in  the  context  of  analysis.  Obviously  there  is  a  great  amount  of  information 
contained  in  the  experimental  results.  The  qualitative  trends  have  been 
discussed  thoroughly  in  Reference  18,  but  it  seems  worthwhile  including 
here  a  brief  summary  of  the  main  conclusions. 

(i)  For  a  given  propellant  and  grain  geometry,  there  is  always  an 
operating  pressure  above  which  finite  disturbances  are  unstable. 
This  has  been  found  true  for  all  propellants  tested  by  both  SRI 
and  GARDE.  It  should  be  noted  that  theoretically  the  effects  of 
the  mean  flow,  and  hence  Mach  number,  are  important  to  the 
question  of  stability.  No  attempt  has  been  made,  apparently 

to  correlate  data  with  this  parameter,  or  to  determine  its  influence. 
The  proper  strategy  is  by  no  means  obvious  at  the  present  time. 

(ii)  Although  the  time-averaged  pressure  increases  with  time 
(see  Figure  4)  the  amplitude  of  the  wave  (change  of  pressure 
divided  by  the  mean  pressure)  measured  at  the  head  end  remains 
essentially  constant.  The  amplitude  is  also  insensitive  to  the  ratio 
of  port-to-throat  area,  and  Thus,  there  appears  to  be  a 

limit  cycle  determined  largely  by  the  processes  occurring  during 
propagation . 

(iii)  Tests  in  motors  geometrically  scaled  from  lengths  of  20  inches 
to  80  inches  show  that  the  larger  motors  remain  nonlinearly 
stable  to  higher  operating  pressures .  The  design  of  the  nozzle 
has  very  little  influence.  Increases  to  larger  sizes  (180  inches) 
showed  only  a  small  further  change  in  the  minimum  stable 
operating  pressure. 

(iv)  For  a  given  propellant  and  grain  cross-section,  motors  having 
"low"  Values  of  length  to  port  diameter  (and  hence  higher 
frequencies  for  the  longitudinal  mode)  tend  to  be  more  stable. 

"Low"  here  means  about  five  or  less.  At  larger  values,  there 
seems  to  have  been  relatively  little  influence  on  the  unstable 
waves.  Whether  or  not  his  behavior  is  closely  related  to  the 
frequency  response  of  the  propellants  is  not  known. 


(v)  Changing  the  grain  cross-section  has  essentially  no  influence 
if  is  fixed.  In  other  words,  the  waves  really  are  close  to 
being  one-dimensional. 

(vi)  Very  roughly,  propellants  having  high  burning  rates  tend  to  have 
higher  threshold  pressures.  Put  another  way,  it  appears  that 

if  all  variables  including  the  mean  pressure  are  fixed,  an 
increase  in  burning  rate  (as,  for  example,  by  adding  rate  modifier 
or  ballistic  additive)  is  a  stabilizing  influence.  This  is  contrary 
to  much  of  the  evidence  for  linear  stability. 

(vii)  Addition  of  lithium  fluoride  which  reduced  the  burning  rate 
of  an  aluminized  propellai.  greatly  reduced  the  threshold 
pressure  and  in  fact  produced  "extremely  violent  instability 
at  practicable  levels"  (Ref.  13). 


(viii)It  is  difficult  to  generalize  about  the  influence  of  aluminum 
(Ref.  18).  However,  it  certainly  is  true  that  the  addition  of 
aluminum  is  by  no  means  a  guaranteed  stabilizing  influence 
and  may  very  well  aggravate  a  problem.  This  has  been  noted 
also  in  other  works  (Refs.  24  and  25,  for  example). 


In  summary,  the  data  reported  in  these  works  is  related  directly  to 
nonlinear  stability.  However,  since  there  is  no  information  concerning 
the  growth  and  decay  rates  of  pulses,  the  measurements  cannot  be  treated 
quantitatively  within  the  analysis  presently  available  or  covered  In  the  work 
reported  here. 


2.2.2  Other  Experimental  Results 


The  early  observations  reported  in  References  25  and  26  contain 
insufficient  detailed  information  to  permit  quantitative  interpretation.  One 
of  the  main  conclusions  of  Reference  25,  however,  is  of  considerable 
practical  interest,  and  is  consistent  with  (viii)  above:  The  addition  of 
aluminum  may  cause  a  stable  motor  to  become  unstable  in  the  longitudinal 
mode. 
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Although  the  observations  discussed  briefly  in  Ref.  24  were  taken  in 
a  center-vented  configuration,  which  will  not  be  treated  in  the  analysis 
discussed  here,  the  data  are  important  for  at  least  two  reasons.  First, 
the  unstable  motions  arise  spontaneously,  and  second  they  exhibit  very 
distinct  nonlinear  behavior. 

The  more  recent  measurements  taken  by  Micheli  (Ref.  23)  are 
potentially  very  useful.  These  comprise  approximately  two  dozen  firings 
of  tubular  motors  having  a  fundamental  frequency  of  about  800  Hz  for  the 
longitudinal  mode.  The  motors  were  pulsed,  but  spontaneous  instabilities 
were  also  observed.  Detailed  records  are  available,  but  unfortunately  the 
data  have  not  been  reduced  from  the  raw  state.  It  is  particularly  Important 
that  growth  and  decay  rates  can  be  found  from  the  records.  The  instabilities 
appear  to  be  mainly  standing  waves. 

The  same  propellant  (ANB3066)  was  used  in  all  tests.  As  a  result  of 
current  efforts  in  several  laboratories ,  the  dynamical  characteristics  should 
be  quite  well  known  in  the  near  future.  One  Interesting  conclusion,  based 
on  a  comparison  of  two  firings,  is  that,  as  Brownlee  has  reported,  an 
increase  of  mean  pressure  set  to  be  a  destabilizing  influence.  Linear 
analysis  (Section  8)  shows  that  .ne  major  reason  for  this  behavior  is  the 
decrease  of  the  Mach  number  of  the  mean  flow,  and  hence  a  decrease  in  the 
losses  at  the  aft  end,  as  the  mean  pressure  is  increased.  Some  of  the 
numerical  examples  presented  later  in  the  present  work  have  been  based 
on  these  pulsed-motor  firings . 

2.3  Analytical  Work  on  Longitudinal  Instabilities 

The  numerical  analysis  discussed  in  the  present  report  is,  in  an 
essential  way,  a  new  contribution  to  understanding  the  problem  of  nonlinear 
instabilities  in  solid  propellant  rocket  motors.  That  is,  in  no  previous  work 
have  the  gas  dynamical  nonlinearities  been  taken  into  account  in  a  computation 
of  transient  wave  motions.  In  this  section,  what  has  previously  been  done 
is  briefly  summarized. 
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Available  analytical  works  which  are  relevant  here  may  be  conveniently 
divided  into  two  classes:  those  which  are  concerned  primarily  with  the  problem 
in  solid  propellant  rockets,  and  those  which  are  concerned  with  other  kinds 
of  nonlinear  wave  motions  mainly  in  liquid  rocket  motors.  None  treat  the 
problem  in  the  detailed  quantitative  manner  discussed  here.  On  the  other 
hand,  they  may  provide  certain  results  more  easily  and  economically  than  an 
elaborate  numerical  calculation - 

A  very  special  feature  of  the  longitudinal  instability  in  solid  propellant 
rocket  motors  is  that  there  are  necessarily  substantial  fluctuations  of  the 
gas  velocity  parallel  to  the  burning  surface.  The  combustion  processes  are 
surely  sensitive  to  these  unsteady  motions,  a  phenomenon  called  "velocity 
coupling."  There  must  also  be  fluctuations  of  pressure,  to  which  the 
combustion  processes  will  respond,  producing  "pressure  coupling." 

Pressure  coupling  must  always  be  present  under  unsteady  conditions,  but 
velocity  coupling  need  not  be,  as,  for  example,  in  purely  radial  modes  in 
a  cylindrical  chamber,  or  in  an  end-burner  oscillating  in  a  longitudinal  mode. 

Consequently,  much  of  the  work  concerned  with  longitudinal 
instabilities  in  solid  propellant  rockets  has  been  directed  to  characteristics 
and  possible  effects  of  velocity  coupling.  (Refs.  27-311..  A  fundamental 
characteristic  of  velocity  coupling  is  that  it  is  intrinsically  nonlinear.  This 
is  a  consequence  of  the  fact  that  the  combustion  processes  respond  to  the 
magnitude  but  not  the  direction  of  velocity  fluctuations,  thereby  introducing 
rectification  effects.  This  idea  has  been  discussed  in  considerable 
detail  in  the  works  referred  to.  The  stability  of  motions  and  the  possible 
generation  of  harmonics  have  been  discussed.  But  growth  of  waves  to  a 
limiting  amplitude,  and  the  nonlinear  effects  within  the  volume  of  the  chamber 
have  not  been  studied.  It  should  be  noted  that  since  the  nonlinearity 
associated  with  velocity  coupling  is  of  first  order  in  the  amplitude,  while 
the  cas  dynamical  nonlinearities  are  of  second  order  (and  higher)  the 
nonlinear  influence  of  the  velocity  coupling  should  be  significant  at 
relatively  lower  amplitudes.  It  is  not  presently  possible  to  establish  this 
conclusion  definitively. 

Two  qualitative  consequences  of  nonlinear  velocity  coupling  are  that 
higher  harmonics  should  be  generated  and  that  an  increase  of  mean  chamber 
pressure  may  be  produced.  Both  features  are  evident  in  some  of  the  data. 
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The  connection  is  appealing  but  cannot  at  this  time  be  supported  quantitatively. 
The  major  problem  in  doing  so  is  that  the  coupling  cannot  be  described 
quantitatively.  Particularly  in  Reference  30,  and  in  subsequent  unpublished 
work,  one  may  find  speculations  and  estimates,  but  all  conclusions  remain 
tentative.  However,  nonlinearities  and  the  mean  flow  in  the  chamber  are 
ignored . 

In  Reference  22,  the  group  at  SRI  attempted  to  interpret  their  data  in 
terms  of  a  model  for  the  coupling  between  a  small-amplitude  pressure 
disturbance  and  the  combustion  processes  at  the  surface.  Certain  aspects 
of  the  gas  dynamics  in  the  chamber  (but  not  the  mean  flow)  are  discussed 
qualitatively  and  to  some  extent  quantitatively.  This  work  was  based  on 
the  idea  that  the  instability  observed  is  a  shock  wave  sustained  by  mass 
addition  at  the  boundary,  which  in  turn  is  related  to  a  response  function 
for  pressure  coupling  only.  They  supposed  that  the  response  function  had 
to  have  a  certain  arbitrarily  chosen  minimum  value  in  order  that  the  wave 
exhibit  sustained  periodic  motion.  They  then  examined  the  values  of 
parameters  appearing  in  their  formulation  of  the  response  function  to 
determine  what  values  were  required  to  match  some  observed  stability  data . 

Although  those  results  appeared  at  the  time  to  be  satisfactory,  they 
are  severely  limited.  No  attempt  was  made  to  account  in  detail  for  the 
losses  in  the  system;  only  a  few  test  results  were  checked.  Moreover, 
the  relationship  of  the  technique,  and  data  obtained  in  the  laboratory  to  the 
full-scale  motors  has  not  been  treated. 

The  problem  of  combustion  instability  in  liquid  rocket  motors  has  of 
course  motivated  a  large  number  of  works  in  the  past  twenty  years.  Of  those. 
References  32-34  are  concerned  with  the  longitudinal,  discrete  wave  motions. 
The  techniques  used  are  applicable  only  to  periodic  shock  waves  and  hence 
have  nothing  to  do  with  the  transient  growth  of  waves,  although  nonlinear 
stability  is  examined.  The  same  sort  of  analysis  was  applied  in  Reference 
35  to  an  end-burning  solid  propellant  rocket  motor.  Obviously,  because  of 
that  restriction,  the  important  influence  of  mass  addition  along  the  lateral 
boundary  is  absent. 

A  modification  of  Galerkin’s  method  has  been  used  in  References  37  and 
38  to  study  stationary  nonlinear  wave  motions  in  liquid  propellant  motors. 


Extension  of  the  technique  to  solid  propellant  motors  with  mass  addition  at 
the  boundary,  and  to  study  travelling  discrete  pulses  is  not  obvious  and 
might  encounter  serious  difficulties.  Moreover,  the  transient  growth  or 
decay  is  not  simply  represented,  and  ultimately  numerical  computations 
are  required.  The  important  advantage  of  the  approach  is  that  only  ordinary 
nonlinear  differential  equations  need  be  treated. 

In  Reference  38,  the  observation  that  the  instabilities  in  solid 
propellant  motors  often  remain  simply  sinusoidal  was  used  as  the  basis  for 
constructing  a  single  nonlinear  ordinary  differential  equation  to  describe 
the  motion.  The  equation  is  that  for  a  nonlinear  oscillator,  and  its  solution 
gives  a  very  good  qualitative — to  a  limited  extent  quantitative— representation 
of  transient  growth  of  waves  in  T-burners.  The  failure  of  that  work  to 
produce  anything  like  the  harmonic  content  actually  observed  (albeit  small) 
has  motivated  the  work  in  Reference  39.  It  appears  that  the  results  will  be 
simple  to  use  and  applicable  to  three-dimensional  as  well  as  one¬ 
dimensional,  motions. 

All  of  the  analyses  discussed  above  involve  approximations  of  one 
sort  or  another.  It  is  practically  impossible  to  f’^termine  how  good  they  are, 
particularly  when  the  periodic  motions  are  not  re.  resented  as  simple  shock 
waves.  Comparisons  with  experimental  data  are  not  entirely  satisfying 
since  usually  fairly  good  results  can  be  obtained  by  choosing  appropriate 
values  for  the  many  parameters  which  necessarily  arise.  The  situation  will 
be  much  improved  if  an  "exact"  numerical  analysis  becomes  available  for 
checking  approximate  analyses.  Thus,  even  if — as  is  likely  the  case —  it 
should  be  impractical  to  use  numerical  calculations  to  attempt  correlations 
of  all  available  data,  nevertheless,  an  important  gap  is  filled.  A 
reasonable  goal  to  head  for  is  an  approximate  analysis  valid  for  any 
configuration,  and  an  accurate  numerical  analysis  which  can  be  used  to 
check  the  approximations  used,  at  least  for  one-dimensional  problems. 

2.4  Nonlinear  Analyses  of  Instabilities  in  Liquid  Rocket  Motors 

In  addition  to  the  works  discussed  above,  there  are  some  calculations 
of  nonlinear  waves  in  liquid  rocket  motors  which  are  not  applicable  to  the 
problem  covered  here  but  should  be  mentioned  for  completeness.  One  of 


the  most  widely  applied  nonlinear  instability  models  was  originated  by 
Priem  and  Guentert  (Ref.  40).  Their  work  has  been  modified  and  extended 
by  several  different  investigators  (Refs.  41-47),  however,  the  model  in  all 
its  forms  is  basically  limited  to  the  conslde»‘ation  of  a  two  phase,  liquid 
drop-gas,  reacting  mixture  in  a  thin  annulus.  Depsite  some  of  the  short¬ 
comings  of  the  model,  it  does  lead  to  the  definition  of  stability  limits  in 
terms  of  engine  design  parameters;  a  feature  which  has  led  to  a  better 
understanding  of  liquid  rocket  instability.  Priem's  model  is,  however, 
not  particularly  well  suited  for  an  extension  to  solid  rocket  engines,  and 
is  limited  to  the  consideration  of  transverse  instability  modes  which  are 
not  being  considered  in  the  current  program. 

Burstein  and  Schechter  (Ref.  48)  have  developed  two-dimensional 
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transient  models.  Two  separate  and" complementary  programs  were  developed 
describing,  respectively,  a  pancake  type  motor  (r-0  model)  and  a  toroid 
with  incremental  thickness  Ar  (toroidal  or  0  -z  model).  The  toroidal  model 
is  the  first  nonlinear  analysis  ior  investigating  both  tangential  and 
longitudinal  motions.  Limited  results  were  obtained  mainly  due  to  the 
excessive  expense  associated  with  solving  two-dimensional  transient 
governing  equations.  The  pressure  amplitude  determined  using  the  pancake 
model  appeared  to  be  unrealistically  high.  This  can  be  partly  attributed  to 
energy  accumulation  in  the  plane  z  =  constant.  This  is  very  similar  to  the 
constant  energy  assumption  in  the  annulus  required  by  the  Priem  analysis. 

Recently  Agosta  (Ref.  49)  developed  a  three-dimensional  transient 
analysis.  Numerical  analysis  was  only  carried  out  only  for  the  one¬ 
dimensional  case.  Nonuniform  droplet  temperature  was  considered  as  a 
result  of  assuming  finite  thermal  conductivity.  Evaporation  kinetics  for 
nonequilibrium  conditions  were  also  included  (Ref.  50)  although  the  validity 
of  this  work  at  high  pressure  has  yet  to  be  proved.  Theoretical  wave  form 
results  were  obtained  for  the  transient  one-dimensional  model. 

Many  nonlinear  extensions  to  the  linear  theory  of  the  time  lag  model 
have  been  developed  in  studying  instability  of  liquid  rocket  engines.  Various 
approximation  methods  were  used  and  the  analyses  do  not  depend  totally  on 
computerized  solutions  as  those  previously  described.  Based  on  this  time-lag 
concept,  Sirignano  and  Crocco  (Ref.  32)  studied  longitudinal  combustion 
instability  for  pressure  waves  of  finite  amplitude.  Mass  and  energy  addition 
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was  assumed  to  occur  only  in  an  arbitrarily  thin  region  next  to  the  injector 
face.  The  solutions  were  shown  to  be  unstable  thus  indicating  the 
possibility  of  triggering  longitudinal  instability.  Based  on  similar  concepts 
Zinn  (Ref.  51)  studied  transverse  mode  instability  using  r>scillatory  nozzle 
flow  conditions  while  Sirignano  employed  the  short  nozzli'  concept. 

Michell  (Ref.  33)  extended  the  work  of  Sirignano  to  include  possible 
discontinuous  waves  and  distributed  combustion.  In  this  manner  he  was 
able  to  show  that  the  final  form  of  triggered  longitudinal  instability 
consisted  of  shock  waves  moving  back  and  forth  along  the  chamber.  As 
noted  above,  the  Galerkin  method,  with  some  modifications  has  been 
applied  to  some  problems  associated  with  combustion  instability  in  liquid 
rocket  engines.  The  important  simplification  achieved  by  the  Galerkin 
method  is  that  the  nonlinear  partial  differential  equations  governing  the 
problem  are  replaced  by  nonlinear  ordinary  differential  equations.  This  is 
achieved  by  direct  integration  of  the  original  equations  over  the  volume  of 
the  chamber.  Time  remains  as  the  sole  independent  variable.  It  is  therefore 
natural  to  consider  a  disturbance  which  is  distributed  in  space  and  maintains 
essentially  the  same  spatial  form  in  time,  but  has  amplitude  which  changes 
nonlinearly  in  time. 

Unfortunately,  the  analyses  developed  for  liquid  rocket  instability 
studies  cannot  be  easily  applied  to  investigate  solid  propellant  motor  problems. 
One  of  the  most  crucial  parts  of  all  combustion  instability  problems  is  the 
coupling  between  the  flow  and  the  combustion  processes.  It  is  there  that 
the  source  of  energy  for  exciting  and  maintaining  combustion  instability 
resides.  Heterogeneous  combustion  processes  of  solid  propellants  are  totally 
different  from  liquid  droplet  burning  which  are  usually  governed  by  the  rate 
of  vaporization  and/or  gas  phase  kinetics,  furthermore,  the  spatial 
distribution  of  solid  propellant  grains,  as  '  .  the  case  of  motors  with 
cylindrically  perforated  prppellant  grains,  presents  an  entirely  different 
problem  formulation  as  compared  with  head-end  injectors  employed  in  liquid 
engines.  As  shown  very  strikingly  by  the  results  obtained  here,  the  truly 
transient  character  of  the  burning  under  unsteady  conditions  is  a  fundamental 
part  of  the  problem.  Essentially  all  of  the  computations  for  instabilities 
in  liquid  rockets  involve  the  assumption  of  quasi-steady  coupling  between 
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the  flow  and  combustion.  If  that  assumption  were  used  in  the  present 
work,  the  results  would  be  unrealistic  and  misleading. 

The  works  of  Powell  and  Zinn  (Refs.  36,37)  are  most  applicable  to 
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the  analysis  of  standing  waves  where  disturbances  are  distributed  in  space 
and  maintain  the  same  spatial  form  in  time.  These  methods,  with  difficulty, 
could  probably  be  modified  and/or  extended  to  treat  solid  rocket  longitudinal 
instability.  However,  it  is  doubtful  that  all  of  the  nonlinear  coupling  effects 
could  be  properly  accounted  for. 
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3.  EQUATIONS  OF  MOTION 

3.1  Discussion  of  the  Mathematical  Model 

The  flow  in  the  combustion  chamber  of  a  metal  loaded  solid  propellant 
rocket  engine  is  calculated  from  a  set  of  one-dimensional,  unsteady,  two- 
phase  flow  equations.  Before  presenting  the  equations,  the  assumptions 
upon  which  they  are  based  are  discussed. 

It  is  assumed  that  the  gas  is  ideal  and  nonreacting;  the  flow  is  inviscid 
and  one-dimensional;  the  particles  are  spheres  of  a  single  size  with 
uniform  Internal  temperature,  the  particles  are  of  negligible  volume  and  do 
not  interact  with  each  other.  The  mass  and  particles  coming  from  the 
burning  surface  are  assumed  to  enter  the  chamber  normal  to  the  burning 
surface  with  zero  tangential  momentum,  and  at  the  adiabatic  flame  tempera¬ 
ture. 

Most  solid  rocket  engine  combrstion  chambers  do  not  have  rapid  area 
changes,  and  the  time  for  a  pressure  signal  from  the  propellant  surface  to 
reach  the  centerline  is  usually  a  small  fraction  of  the  period  of  the  longi¬ 
tudinal  pressure  waves;  hence,  a  one -dimensional  analysis  should  provide 
a  reasonable  approximation  to  the  flow. 

The  combustion  of  a  metal  loaded  solid  propellant  does  not  directly 
produce  solid  particles.  The  solid  particles — aluminum  oxide  for  aluminized 
propellants — are  formed  by  the  combustion  of  metal  droplets  which  have 
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been  entrained  by  the  gas  flowing  over  the  burning  surface'  .  In  the  current 
work,  however,  the  formation  and  combustion  of  the  metal  droplets  has  not 
been  modeled  and  inert  solid  particles  are  assumed  to  be  carried  into  the 
flow  by  the  burning  gases  at  a  specified,  and  constant  weight  fraction. 

In  the  analysis  the  particles  are  also  assumed  to  be  all  of  the  same 
diameter,  and  all  particles  at  a  given  location  are  assumed  to  have  the 
same  velocity  and  temperature.  In  actuality,  the  particle  sizes  form  a 
distribution,  which  is  at  least  bimodal  due  to  two  modes  of  oxide  for¬ 
mation  Quantitative  particle  distribution  functions  are  not  available, 
but  the  particle  formation  mechanisms  are  such  as  to  produce  both  small 
0-2^particles  (smoke)  and  larger  5-20*1  particles. 


Particle  injection  takes  place  continuously  along  the  entire  length  of 
the  chamber.  Within  a  small  volume,  at  any  one  time,  various  particles 
may  possess  different  velocities  and  temperatures  (even  if  they  are  all  of 
the  same  size).  This  is  a  result  of  the  fact  that  the  history  (trajectory)  of 
each  particle  is  different,  and,  in  time,  particles  from  different  locations, 
injected  at  different  times,  can  be  in  the  same  place.  The  ability  to 
distinguish  between  particles  because  of  their  varying  past  histories  is 
diminished  with  time,  since  all  the  particles  eventually  acquire  the  local 
velocity  and  temperature.  For  particles  of  radius,  a  ,  mass,  m,  moving  in 
a  gas  with  viscosity,  u  »  characteristic  times  required  for  the  particle  to 
approach  the  local  velocity  and  temperature  may  be  defined  (for  Stokes 
flow) 
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Under  typical  conditions  encountered  in  a  solid  rocket  motor  and 
are  usually  about  0.1  to  0.01  milli-seconds;  only  a  fraction  of  the 
typical  period  of  a  longitudinal  wave.  Therefore,  the  current  assumption 
that  all  particles  at  a  given  location  and  time  have  the  same  velocity 
and  temperature,  is  a  reasonable  one.  Although  this  assumption  is  a 
reasonable  one  it  raises  the  conceptual  problem  of  how  the  particles 
injected  at  a  given  time  instantaneously  acquire  the  local  velocity  and 
temperature  of  the  gas.  There  are  two  modes  by  which  this  nonphysical 
event  could  be  postulated  to  occur.  The  required  energy  and  momentum 
could  be  instantly  acquired  from  the  gas,  corresponding  to  infinite  drag 
coefficient  and  Nusselt  number.  Or,  alternatively,  the  newly  injected 
particles  may  be  said  to  instantly  exchange  energy  and  momentum  with  the 
other  particles  through  collisional  effects  (infinite  cross-section) .  The 
choice  of  one,  or  the  other,  of  these  postulates  is  reflected  in  the  momentum 
and  energy  equations;  however,  as  might  be  expected,  as  long  as  r^and 
T.J.  are  small  compared  to  typical  particle  stay  times  in  the  chamber,  the 
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choice  makes  little  difference. 

The  equations,  as  used  herein,  imply  the  assumption  of  instantaneous 
acquisition  of  momentum  and  energy  through  collisions.  The  alternate  forms 
of  the  equations  are  presented,  for  reference,  in  a  footnote  to  the  equations. 

The  assumptions  previously  mentioned,  and  discussed  were  made  to 
allow  the  formulation  of  a  relatively  realistic  nonlinear  combustion 
instability  model,  without  undue  complication.  The  present  model  should 
allow  for  a  reasonable  assessment  of  this  type  of  approach,  and  should  it 
prove  to  be  warranted,  the  model  may  be  enhanced  by  the  removal  of  one, 
or  all,  of  the  above  assumptions. 

3.2  Conservation  Equations 

Subject  to  the  previous  discussion,  the  equations  of  motion  for  the 
two  phase  mixture  are  as  follows: 
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where  asterisks  denote  dimensional  variables,  p*  and  p*  are  the  gas  and 
particle  densities,  respectively,  u*  and  u*  are,  respectively,  the  gas 
and  particle  velocities.  A*  is  the  chamber  cross-sectional  area,  while 
uu*  and  uu*  are  the  mass  flux  of  gas  and  particles  given  off  by  the  propellant 
per  unit  length,  per  unit  cross-sectional  area. 

The  quantity,  u)*,  is  related  to  the  burning  rate  of  the  solid  propellant, 
r*,  as  follows; 


P*r*je* 

,*  =  ■?-TT- — 


A*(l+p^) 


(3-4) 


i**'*i*t*. 


where  p*  is  the  density  of  the  solid  propellant,  X*  is  the  perimeter  of  the 
s 

grain,  and  is  the  weight  ratio  of  solid  particles  to  gas.  The  steady 
burning  rate,  r*.  Is  assumed  to  have  the  form 


Momentum 

Gas  Phase  ^  (p*u*A*)  +:^-*(p*u*^A*)  =  -  mA*+  F*A*  (3'«) 
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Particles  ^  (pJuJA*)  (PjuJ^A*)  =  -  F*k*  (3-7) 


where  the  term  F*  represents  the  effect  of  momentum  transfer  between  the 
particle  and  gas.  The  momentum  interaction  is  equal  to  the  drag  force 
exerted  by  the  particles  on  the  gas,  per  unit  voiume,  and  is  given  by 


r,  9*9* 
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where  a*  is  the  particle  radius,  p*  the  density  of  the  solid  particles  and 
Cjj  is  the  drag  coefficient  of  the  particles.  The  drag  coefficient  Cjj  is 
obtained  from  a  well-established  correlation  of  the  so-called  "standard" 
drag  coefficient  versus  Reynolds  number  data .  These  data  represent 
measurements  for  a  single  sphere  in  steady  flow.  The  formula  used  was 

/  go\ 

originated  by  Kliachko'  ,  and  is. 


(3-9) 


The  second  term  in  the  above  equation  represents  a  correction  to  the  Stokes 
value,  24/Re,  and  allows  the  formula  to  be  used  at  Reynolds  numbers  up 
to  several  hundred. 


The  gas  and  particle  momentum  equations  (3-6)  and  (3-7)  can  be 
rewritten  in  modified  form  by  subtracting  the  respective  continuity 
equation,  (3-2)  or  (3-3),  after  it  has  been  multiplied  by  velocity  (u*  or 
u*  as  the  case  may  be).  The  resulting  modified  momentum  equations  are: 


Gas  Phase 


p*  Au?  +  p*u*  ^ 

P  bt*  ^  ^  dx* 


l*  (U  *  — 


Particles 


p  at 


*  a  u* 

§  +  P*U*  ^ — I 

*  ^p  p  a3c* 


u*  tt)*  -  F* 
P  P  P 
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Energy 

The  energy  equations  for  the  particles  and  gas  can  be  written  in  several 
different  ways,  in  terms  of  internal  energy,  enthalpy,  temperature,  entropy, 
etc.  The  equations  are  written  below  in  terms  of  temperature. 


+As  discussed  in  Section  3,1,  the  momentum  and  energy  equations  are  slightly 
altered  if  it  is  assumed  that  the  entering  particles  immediately  acquire  their 
momentum  and  energy  from  the  gas  flow.  The  differences  are  as  follows:  the 
tern,  u*  of  must  be  subtracted  from  the  r.h.s,  of  Eq,  (3-10),  and  added  to  the 
r,h.s,^of%q.  (3-l^)^‘  the  terms  u*  u*  and  w^C|(T^-l^)  -  1/2 u^J 

(tf-  C*{Tt  -  T*)  +  — ^  Jmust  be  subtracted  from  the  r.h.s.  of  Eq.  (3-13). 

p—  S  r  p  4 
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v/here  C*  and  C*  are  the  heat  capacities  of  the  gas  and  solid, 

P  ® 

respectively,  and  Q*  represents  the  volui.  etric  rate  of  heat  transfer 
between  the  particles  and  gas. 


p*  G*  (X* 

P  P  P 


m*C* 

_ 

2NuTTo*kV 


PJC*(T*-T*) 


(3-14) 


where  k*  is  the  gas  thermal  conductivity,  Nu  is  Nusselt  number  and  m*  is 
the  mass  of  a  particle. 


1*  =  1  rrc*^  p* 
3  m 


(3-15) 


The  parameter  =  (m*C*)/(2Nuno*lc*)  is  the  thermal  relaxation  time 
constant,  i.e. ,  is  the  characteristic  time  it  takes  the  particles  to  adjust 
to  the  local  gas  temperature.  The  Nusselt  number  is  a  function  of  Reynold r 
number  and  Prandtl  number.  The  following  expression  suggested  by 
Carlson'®^^  is  used: 

Nu  =  2  +  .459  Re’^^Pr'^^  (3-16) 


Here  again,  the  second  term  represents  a  correction  to  the  Stokes  flow 
value,  i.e.,  Nu  =  2.  The  gaseous  equation  of  state  is  also  required  and  is 

P*=c*RT*  (3-17) 


Nondimensional  Equations 


In  an  effort  to  improve  the  accuracy  and  reduce  the  possibility  of  round¬ 
off  error  in  the  solution,  the  independent  and  dependent  variables  have  been 
normalized.  The  resulting  normalized  (nondiniensional)  variables  are  more 
uniform  in  magnitude.  The  following  new  variables  have,  therefore,  been 
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defined: 


=  xVL 

T  =  TVT* 

A  =  A*/L*^ 

=  t*a*p/L* 

P  =  P*/P|! 

ou  =  U)*L*/p|a* 

=  u*/a* 

p  =  ()*/p* 

=  t«^I-*/p|a*  ^3 

=  pj/pj 

T„  =  T*/T* 

P  P  F 

Fp=  FJLVpfa*2 

=  uVai 

P  F 

C  =  c*/c* 

s  p 

Op=:  OJL*/pJCJa*T 

=  a*/a* 

where  is  the  propellant  adiabatic  flame  temperature,  is  chamber 
pressure,  and  a^  are  the  density  and  sound  speed  evaluated  at  the  chamber 
conditions  and  L*  is  a  reference  length,  usually  the  length  of  the  grain. 

In  terms  of  these  nondimensional  variables  equations  (3-2),  (3-3), 

(3-10  to  13)  and  (3-17)  become: 

Continuity 


^  (PA)  +  (puA)  =  u)A 


(3-i9) 

(3-20) 


Momentum 

~  5  u  1  ^  .  r> 

+  P  u  T”  =  -  UtlJ  “  —  T”  +  F 
^dt  ^  p 


3u  3u 

n  - E  +  D  U  - ^  =  -  («  u  -  F 

3t  Pp^p  ax  p“p  u 


(3-21) 

(3-22) 
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Energy 
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The  momentum  equation  (3-21)  is  then  multiplied  by  an  unknown 
multiplier,  q,  and  added  to  equation  (3-27),  this  gives 


^  ^  Mj  ^  If  ^  ^  -v-^  If 


(3-28) 


"  lx  Fpl_(Y-  l)(u  -u)  +  qj  +  0 


If  the  total  derivatives  of  pressure  and  velocity  are  to  be  the  same 
then  the  condition 


(3-29) 


must  be  satisfied.  Since  a'^  =  P/p  ,  this  condition  leads  to 


q  =  +  a 


(3-30) 


"  ft "  It  ^ 


It  t; 


are  defined  as  the  total  derivatives  along  the  characteristics  lines  given  by 


f  =u±a 


(3-32) 


then  equation  (3-28)  with  (3-30)  yields  the  following  compatlbiity  relations: 


1  6'*’P 

Yoa  6t 


X  S  u  _  a  f  1  3A  .  SA]  .  tt'  T/,  \  ’ 
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oa 
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The  gas  energy  equation  (3-23)  is  already  in  characteristic  form. 
The  gas  streamlines 


(3-35) 


are  characteristic  lines,  and  if 

_ ^ 

6t  ^  5x 


(3-36) 


is  defines  as  the  total  derivative  along  the  streamlines ,  then 

the  energy  equation  becomes  the  compatibility  relation,  and  may  be  written 

as 


0 


6J 
6  t 


(>  -1)  IP  ^  (v-1)  P  iA  . 
'  Y  A  c' t 


oui'd-T)  + 


+  (v  -  1)  F  (u  -u)  +  Q 

p  p  p 
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A  particle  path. 


dx 

dt 


(3-38) 


is  a  dual  characteristic.  Defining, 


+  u 

P 


9  X 


(3-39) 
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as  the  total  derivative  along  a  particle  path,  allows  the  particle  momentum 
and  energy  equations  (3-22)  and  (3-24)  to  be  written  as 


0 

P 


6  u 

-E_P  = 
6t 


P  P  P 


(3-40) 


6  T 
-E_a  = 


^p  6t 


=  ou  (1  -T  ) 
P  P 


/v-1 

V  2 
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The  particle  continuity  equation  (3-20)  becomes  uncoupled  from  the 
remaining  equations  (since  there  is  no  particle  equation  of  state  which 
relates  particle  dansity  to  particle  temperature)  and  cannot  be  written  in 
characteristic  form.  It  can,  however,  be  written  in  the  following  alternative 
form 


5  0 


■J-P 

6t 


^p 


_ P  _E_ 

A  fit 


(3-42) 


The  required  equations  of  motion  are,  therefore:  the  equations  of  the  four 
characteristic  lines,  equations  (3-32),  (3-35),  (3-38);  the  five 
compatibility  relations,  equations  (3-33),  (3-34),  (3-37),  (3-40),  (3-41); 
and  the  particle  continuity  equation  (3-42).  These  equations,  together  with 
the  equation  of  state,  (3-25),  the  transient  burning  rate  analysis  (Section  4), 
boundary  conditions  (Section  5) ,  and  initial  conditions  (Section  6) ,  form  the 
complete  mathematical  model  of  the  instability  problem. 


/l.o  RESPONSE  FUNCTION  FOR  TRANSIENT  BURNING 


In  general  one  must  expect  that  the  combustion  processes  will  respond 
to  both  pressure  and  velocity  fluctuations.  Moreover,  although  this  has 
not  been  established  experimentally,  the  response  may  be  nonlinear  in 
addition  to  the  Intrinsic  dependence  on  the  magnitude  of  velocity  fluctuations 
mentioned  in  Section  2.3.  Other  than  a  numerical  calculation  of  tlie 
response  to  an  exponential  change  of  pressure  (Ref.  58),  there  have  been 
three  treatments  of  the  nonlinear  response  to  pressure  coupling  (Refs,  55,56, 
57),  All  are  restricted  to  harmonic  motions,  and  seek  to  obtain  nonlinear 
heat  conduction  solutions.  There  are  open  questions,  however,  about  the 
validity  of  these  analyses.  The  present  treatment  Is  based  on  a  linear 
solution  for  the  heat  conduction  in  the  solid,  coupled  to  a  nonlinear  solution 
of  the  chamber  flow  dynamics.  This  approach  seems  to  be  justified  based  on 
observation  of  data. 

The  method  for  calculating  both  the  pressure  coupled  (Section  4,2)  and 
velocity  coupled  (Section  4.3)  response  to  disturbances  of  arbitrary  waveform 
is  essentially  an  extension  of  the  analysis  of  the  response  to  harmonic 
oscillations.  Therefore,  a  summary  of  the  harmonic  analysis  is  included  here 
to  clarify  the  origin  of  the  parameters  required  in  the  numerical  analysis.  The 
influence  of  aluminum  is  not  explicitly  accounted  for, 

4,1  Linear  Response  To  Harmonic  Pressure  Oscillations 

It  has  been  shown  (Ref.  59)  that  almost  all  existing  analyses  of  the 
response  of  a  burning  solid  to  harmonic  pressure  fluctuations  lead  to 
essentially  the  same  functional  dependence  on  frequency.  The  reason, 
logically  enough,  is  that  all  are  based  on  the  same  assumptions,  namely: 

(1)  The  gas  and  solid  phases  are  treated  as  homogeneous 
phases,  and  only  variations  in  the  direction  normal  to  the 
burning  surface  are  accoi  nted  for. 

(2)  Conversion  of  solid  to  gas  occurs  at  an  infinitesimally  thin 
interface;  in  particular,  no  chemical  reactions  are  present 
in  the  solid. 

(3)  The  gas  phase  responds  quasi-statically  so  that  only  the 
steady-state  solution  is  required.  This  assumption  is 
justified  by  the  fact  that  processes  in  the  gas  phase 

are  much  faster  than  those  in  the  solid. 
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The  analysis  is  then  conveniently  split  into  three  parts:  the 
solid,  the  gas,  and  the  interface. 


4.1.1  Solid  Phase 

A  coordinate  system  fixed  to  the  gas/solid  interface  is  used.  With 
the  approximations  listed  above,  and  assuming  also  that  the  material 
properties  are  constant,  the  problem  reduces  to  the  heat  conduction 
equation, 


D  r  ^ 
^s^s?  t 


OoC 

s  S 


11  ^ 

?x 


2 

IIT 


(4 


dx' 


The  equation  is  linearized  by  writing  T  =  T  +  T* ,  r  =  r  +  r';  the  steady- 
state  and  fluctuating  temperatures  satisfy 


s  dx  "  ^s  .  2 
dx 


p  C  P-'  +  p  CF  ^’=k  ^ 

s  s8t  ^s  /  ;^x  .  2 

c  X 
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From  (4-2),  with  the  conditions 


T  =  Tg  (x  -  - 


T  =  T  (x  =  0) 

w  '  ' 


(4 


(4 


(4' 


one  finds 


It  is  sufficient  for  linear  problems  to  examine  the  case  of  harmonic  time 
variations.  The  solution  to  Eq.  (4-3),  yields  the  formula  for  fluctuations 
of  heat  transfer  from  the  interface  into  the  solid, 


=  mCj  XT’ 
s.  w 


(4-7) 


where  m^  =  p^r'  is  the  fluctuation  of  mass  flux,  T|^  is  the  fluctuation  of 
surface  temperature,  and  X  is  the  complex  function  of  frequency: 


X  =1/2  'j.  +  ^  \  -/T+iTn^  +  1  j 


(4-8) 


H 

n  =  -^ 


U) 


(4-9) 


where  uj  is  angular  frequency  of  the  waves. 

It  is  common  practice  to  assume  that  a  simple  pyrolysis  law  of  the 
Arrhenius  form  applies: 


n  -E  /R  T 
=  B  p'^e  w  ow 
w  w^ 


(4-10) 


Thus,  the  fluctuations  of  mass  flux  and  surface  temperature  are  related  by 


S-  =E  ^  +  n 
m  =  w  p 

^w 


(4-11) 


with  the  normalized  activation  energy  for  the  surface  reaction 


E  =  E  /R  T 
wow 


(4-12) 


Equations  (4-7)  and  (4- 11) may  be  combined  to  yield  the  desired  result. 


L  s  a  xj  - 


(4-13) 
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where  the  important  parameter  A  is  defined  as; 


A=^;i--=)E 

w 


(4-14) 


4.1.2  Interface 

By  examining  conservation  of  mass  and  energy  for  a  control  volume 
encompassing  the  interface,  one  finds  the  general  "jump"  conditions 


s  ►^g  g 
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in  which  the  enthalpies  on  the  solid  and  gas  sides  of  the  interface  are 
h  , ,  respectively.  The  heat  of  reaction  associated  with  the  processes  at 

W 


the  interface  is  Q,.,, 

W 


^w  ^w-  ”  ^w+ 
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and  Q  >  0  for  an  exothermic  reaction, 
w 


The  perturbation  of  (4-16)  is 


HlV'h  fi]'-  ""■pp-Cs>V'”w'?v 


(4-18) 


and  with  (4-13),  one  has 


C  T' 
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where 


H  = 


°s<fw-V 


(4-20) 


4.1.3  Gas  Phase 

This  is  the  most  difficult  part  of  the  problem  if  one  wishes  to  cope 
with  the  details  of  the  flame.  However,  it  turns  out  that,  owing  to  the 
assumption  of  quasi-static  behavior,  no  matter  what  model  one  selects  for 
the  flame  structure,  the  ultimate  form  of  the  response  functions  will  be  the 
same.  For  simplicity,  the  most  elementary  model  will  be  discussed. 

The  energy  equation  for  a  thermal  theory  of  flame  propagation  is 


mCp 


dx  ~  dx  vg  dx>‘ 


+  Q{  w 


(4-21) 


where  Q.  is  the  heat  release  per  unit  mass,  and  w  is  the  reaction  rate 
-1  ^ 

(sec  )  divided  by  the  gas  density.  This  equation  must  be  solved  subject 
to  the  conditions 


(4-22) 


at  the  surface  (x  -•  0) ,  and 


T  =T 


f 


dx 


=  0 


(4-23) 


at  the  downstream  edge  of  the  flame  (x  -♦  <»). 

Now  suppose  that  the  energy  release  is  uniform  in  the  region  from  the 

surface  to  the  edge  of  the  flame  (x  =  x^)  so  Q^w  is  constant.  This  case  has 

been  treated  in  References  58  and  59.  Equation  (4-21)  can  then  be  integrated 

directly  (note  that  m  =  p  u  is  constant)  to  give  the  formula  for  the  heat 

9  ^ 
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transfer  to  the  surface 
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It  happens  that  for  typical  values  of  the  quantities  involved,  the  exponential 
involving  is  negligible: 

mC  p  r  _  C  k  .. 

- EL  X  =  C  X  =  I  --B-S.  1-^  X 

^  P  f  Lc^kg  J  ^f 
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(1  cm/sec)  ^ 
(10~^  cm/sec)  ^ 


1000  x^ 


Even  if  Xj  is  as  unrealistically  small  as  20  microns,  1000  x^  =  2;  hence  the 
approximation  is  quite  good,  and  (4-2^  becomes 


Ik  =  j'£f-g'|  w 

LgdxJ+  L  C  Jm 


The  pertu’^ba tion  of  the  last  equation  is 
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It  remains  only  to  find  an  expression  for  w,  and,  hence,  w', 
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(4-25) 


(4-26) 


(4-27) 
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^  Suppose  that  the  reactions  are  mainly  sensitive  to  changes  of  pressure, 
so  w'  ^  p* .  It  is  convenient  to  define  Why 
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s  w 


With  this  assumption,  (4-26)  becomes 


•J»  X  p 
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(4-29) 


which  is  the  result  required  of  the  solution  for  the  gas  pha 


se. 


4.1.4  Linear  Response  Function  and  Fluctuations  of  the  Flame 
Temperature 


Equations  (4-11),  (4-19),  and  (4-29)can  be  combined  to  yield  the  ratio 

C. 
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pVF~ 
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(4-30) 


In  the  limit  of  zero  frequency,  X  -»  1  and  (4-29)must  reproduce  the 
small  change  of  burning  rate  due  to  a  small  change  of  pressure  in  steady 
state  burning.  For  the  common  case  r  ^  p"^,  this  implies  that  the  ratio 
(4-30)must  equal  n  for  u;  -♦  0;  i.e. , 


••  c  c  c 
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Define  the  parameter  B  as 


(4-31) 


B  = — r  w + ^  1 


(4-32) 


4-7 


and  the  condition  (4-31)  is 


EA^-IiA  +  -^  -  1  +  (A  +  1)  =  AB 
s 

Hence,  (4-30)  has  the  simpler  form 

,  nAB  +  n  (X  -  l) 

m'/m  _  _ w  ^ 

\+  j  -  (1  +A)  +AB 


(4-33) 


(4-34) 


(Incidentally,  with  a  bit  of  effort,  one  can  show  that  (4-33)  implies 
satisfaction  of  the  mean  energy  balance  at  the  gas/solid  interface  written 
in  a  slightly  obscure  form) . 

Equation  (4-34)  is  proportional  to  the  response  function,  Rj^,  usually 
defined  as 
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A  second  quantity,  the  admittance  function,  is  also  important  in  studies 
of  instability.  This  is  defined  as 
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Now  from  the  definition  of  m  =  pu* 
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(4-37) 


*Note  that  owing  to  the  way  in  which  the  mass  flux  is  used  here,  the  assump 
tion  is  implicit  tha^the  propellant  contains  no  metal.  When  applied  to  metal 
lized  propellants,  m  and  m'  must  then  be  the  values  for  the  gases  only. 


With  the  perfect  gas  law,  (4-37)  can  be  written 


T'/T 

Jiere  T'  is  really  T^,  which  can  be  written  as  the  sum  of  an  isentropic 
part  and  a  nonisentropic  part: 


Tf  ^Tj^lsent. 


Thus,  (4-38)  becomes 
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AT*/  Tj 
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The  last  result  is  generally  true  for  linear  variations.  With  the 
analysis  leading  to  (4-34)  one  can  deduce  an  explicit  expression  for  aT^  and 
hence  for  the  admittance  function.  Consider  first  a  control  volume  extending 
from  the  solid/gas  interface  downstream  beyond  the  combustion  zone;  the 
energy  balance  for  this  region  is 


>  =  ""tSf-CpITf- V] 


(4-41) 


where  is  again  the  energy  released,  per  unit  mass  of  flow,  by  chemical 
reactions.  The  perturbu:‘on  of  (4-41)  is 
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The  mean  energy  balance  for  a  control  volume  extending  from  deep  within 
the  solid  to  a  plane  on  the  gas  side  of  the  interface  is 


[''9SJ-h'=s«w-’’s)-Qw 

With  (4-41),  one  has  the  identity  for  use  in  (4-42): 


(4-43) 


Qr-C  (T,-T)  =  C(T  -T)-Q 
f  p  f  w'  s'  W  S  V7 


(4-44) 


Also,  combination  of  (4-25),  written  for  steady  burning,  and  (4-43)  gives 
a  formula  for  : 


^EA^  =  A(1  -H) 


(4-45) 


Now  after  the  Ai*rhenius  law  (4-11)  is  used  to  eliminate  T'  from  (4-42), 

s 

and  Equations  (4-43)  and  (4-45)  are  substituted  appropriately,  the  result  is 


If  If  P  P 


The  bracketed  terms  multiplying  m’/ m  can  be  rewritten,  using  (4-33)  and 
(4-45)  as 


p  “  p 


(4-47) 


and  with  (4-35),  one  now  has 


T*  C  T 

-I)  ^ 

m  w  m  VH  D  P 

Tj  p  Tj 


(4-48) 


To  find  the  nontsentropic  part  of  (4-48),  write  Tf/Tf  as  in  (4-39),  and 
the  last  equation  produces 


T. 


C  =  E  Vvn'S  \  I  V  ^  J  P 


(4-49) 
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This  completes  the  elementary  linear  analysis  for  harmonic  motions. 

The  fluctuations  m',  T^,  etc. ,  are  of  course  complex  amplitudes  with  phases 
measured  with  respect  to  the  pressure  oscillation.  It  should  be  noted  that 
(4-50)  is  a  general  result,  whereas  (4-49)  involves  the  assumptions  listed  at 
the  beginning  of  this  discussion. 

Although  the  flame  temperature  must  in  principle  oscillate  according 
to  the  preceding  result,  this  behavior  is  not  currently  included  in  the  numerical 
analysis.  There  is  no  reason  why  it  cannot  be  incorporated,  and,  in  fact, 
variable  flame  temperature  will  be  added  to  .the  model  as  soon  as  possible. 

In  the  meantime,  the  exclusion  of  this  effect  should  not  affect  the  immediate 
task  of  assessing  the  possible  value  of  the  current  approach  to  solid  rocket 
instability. 


The  response  function.  Equation  (4-34),  has  been  derived  as  the 

sensitivity  to  harmonic  oscillations.  However,  simply  by  taking  iO  as  the 

variable,  (4-34)  can  be  treated  as  a  Laplace  transform.  Henceforth,  the 

influence  of  surface  reactions  will  be  ignored  (n„,  =  0)  so  with  s  =  i  Q , 

w 


nAB  m 


Xls)  +  (AB-A-i)  X(s) 


(4-51) 


Note  that  s  =  iniu  /r*  is  the  dimensionless  Laplace  transform  variable  as 
associated  with  the  dimensionless  real  time  t , 


(4-52) 


4-11 


.  .'■  . . -.  ^ 

. . . .  I 


Thus,  for  any  function  f(t )  having  transform  F(s), 


F(s)  =  ^f(t)  e‘®^  dr 
0 

f(r)  =  ^  ^F(s)e®^  ds 


Now  define  a  new  Laplace  transform  variable  p. 


p  =  1  +  4s 


so  that  the  quantity  X  is 


\  =  I  [l  +  a/1  +4s  j  =  -I  *1  +'fp  ] 


(4-53) 


(4-54) 


(4-55) 


(4-56) 


Substitute  (4-55)  into  the  Laplace  transform  pair  (4-53)  and  (4-54)  to  find 


F(p)  =  \  ^4e"f(T)  dr 


[4e^£(T)J  =  ^  F(p)  ^Pt  dp 


(4-57) 


(4-58) 


where 


-  =  t  /4  =  ~ — * 
■  ^  4h 


(4-59) 


Hence,  what  appears  now  is  4e^  times  the  desired  function  of  time,  not  the 
function  itselt. 

The  inverse  transform  of  (4-51)  gives  m'/  m  as  a  function  of  t: 


1  ^ 
nAB  m 


^  /ri  =  -L.  C- 

m  ’  2tt1  \  - 


>d± 


+{AB-A-1)X+AJ  P 


^(s)  e®  ^  ds  (4-60) 
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and  if  the  substitutions  t  -♦  t,  s  ->  p  are  made, 


~J—  r4eT  ^/  \  i  ^  '• - -  iL 

2nftB  L'*®  m  i,VF-Ci)(V5-<,2)  ^ 


1  r  (!+'> 


■^e^^dp  (4-61) 


_^l 


where  aj  <  02  are  the  roots  of  the  denominator, 


aj  =  -  A(B  -  1)  +  i  \/4A  -(AB  -A  -  1)® 

Oi  =  ^2  =  -  A(B-1)  -  i74A-(AB-A-l)* 


{4-62a) 


(4-62b) 


Thus,  (jj/  02  are  complex  conjugates. 

% 

In  principle,  (4-61)  can  be  used  for  numerical  calculations.  For  by 
use  of  the  convolution  theorem,  if  p'/  p(t)  is  the  inverse  transform  of 
p'/  p(p)  and  K(t)  is  the  inverse  of  the  remainder  of  the  integrand,  then 
formally 

Sjkr"®’  r<-']  =  I  4e?K(e)  d5 


8ife  ■  S  KCj)  |^(t-5)  d; 


(4-63) 


Note  that  the  factor  4e'  must  be  carried  along  according  to  earlier  remarks. 

Unfortunately,  the  function  K(§)  turns  out  to  be  extremely  awkward 
to  handle  numerically.  It  contains  terms  involving  complex  error  functions 

n®j(T-§)  - - 

and  exponentials  such  as  e  erfc(ojV  T-f )  .  This  combination  can 

be  shown  to  be  equivalent  to  the  W(z)  function  defined  on  page  297of 
Reference  70 . 
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Several  different  series  representations  may  be  used  to  represent 
the  above  term.  Each  series,  however,  is  either  convergent  in  only  part 
of  the  complex  plane,  or  if  convergent  everywhere,  is  only  very  slowly 
convergent  in  part  of  the  range  of  interest.  It  turns  out,  therefore,  that 
even  by  using  different  series  with  overlapping  regions  of  convergence 
there  are  parts  of  the  complex  plane  of  interest  where  many  terms  have  to 
be  taken  in  the  series . 

Since,  in  an  instability  solution,  the  turning  rate  integral  must  be 
carried  out  from  0  to  the  current  time  at  each  of  many  x  locations,  every  time; 
the  computation  time  involved  in  evaluating  the  integrand  becomes  prohibitive 
when  the  series  are  so  slowly  convergent.  The  fact  that  complex  arithmetic 
must  be  used  further  serves  to  amplify  the  problem. 

The  numerical  difficulties  of  the  direct  approach  can  be  avoided  by 
use  of  an  alternate  procedure.  If  (4-51)  Is  written  in  terms  of  P,  using 
(4-56)  and  (4-62),  equation  (4-51)  oen  be  written  in  terms  of  P  as, 


|l(p)  =  j'i  +a_-: 


(4-64) 


Now  take  the  inverse  transform  to  find 


Oi+Oo  J,  , 

-  i  Z  (  m 


^  ^  r  HT/  \  -5  -T  r  m*  e 

^  J  m  \/7- ■  \ 


T 

I  1  r 


(4-65) 


Although  this  is  not  an  explicit  formula  for  m'/  m,  the  functions  under  the 
integrals  are  easy  to  handle  numerically. 

Equation  (4-65)  does  present  a  problem  at  the  lower  limit  of  the  integral, 
?  -  0,  sines  the  integrand  becomes  infinite.  This  singuiarity  can  be  handled  ’ 
by  dividing  the  integral  into  two  pieces,  and  approximating  the  part  near  n 


Thus, 


(4-66) 


where  6  0.  If  only  the  first  two  terms  in  a  Taylor's  series  expansion  are 

retained,  the  first  integral  has  the  following  value: 


6  6  6 


(4-67) 


where  powers  of  6  greater  than  3/2  have  been  ignored. 

All  the  Integrals  in  (4-6^  can  be  treated  in  this’ way,  and  after  some 
rearrangement,  the  final  result  used  in  the  numerical  calculation  is 

(t)  [l  -  I  +  01  Ojs] 

'*•  V  TT 


t-6 

s 

0 


=  \  |^(f)  [7=^= 


„  I  “(t-?) 
-  O1O2 j  e 


de 


§) 


T  "  6 


2nAB  [  ^(F)r-7===-  +  d? 

^  P  L  ^/TT(T-5) 


(4-68) 


13/2  £L(,)r2^r2^'T-f  6'/2)^2nAB6 

,  77  3  drVmyp^TT-  3  ^ 


The  accuracy  of  this  formula  has  been  checked  by  examining  several  special 
cases  which  can  be  worked  out  exactly  (see  Appendix  A). 
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The  result  (4- 68) contains,  in  addition  to  the  pressure  index  n,  two 

parameters,  A  and  B.  0\/er  rather  broad  ranges,  they  can  be  chosen 

independently.  Previous  experience  has  shown  that  A  is  probably  in  the 

range  5  <  A<50  while  B<  1.  This  constraint  on  B  assures  (see  Ref.  60  ) 

that  heat  conduction  to  the  solid  from  the  gas  is  positive.  However,  A  and 

B  cannot  be  set  entirely  independently  for  the  following  reason.  Both  (4-61) 

and  (4-68)  produce  a  response,  m'/  m,  which  consists  of  transient  oscillations 

and  a  long-time  "steady-state"  behavior.  It  is  necessary,  for  acceptable 

solutions,  that  the  transients  decay.  This  condition  is  met  if  and  only  if 

(59) 

A  and  B  satisfy  the  inequality  . 

A  <  (4-69) 

One  would  like  to  choose  A  and  B  to  match  the  dynamical  behavior 
of  the  particular  propellant  used.  In  principle,  data  taken  in  T-burners, 
for  example,  should  provide  adequate  basis.  However,  the  available 
experimental  results  are  limited.  Hence  one  is  forced  to  make  crude 
estimates.  One  aid  in  this  respect  is  the  result  that  the  response 
to  harmonic  oscillations  exhibits  a  peak  occurring  approximately  at  the  value 
of  dimensionless  frequency  ,  Tz , 


(4-70) 


4.3  Velocity  Coupled  Response 

At  the  present  time,  the  probelm  of  velocity  coupling  is  unsolved. 

There  is  little  question  that  the  burning  of  a  propellant  will  respond  to 
disturbances  of  the  velocity  parallel  to  the  surface.  For  steady-state 
burning  this  is  called  erosivlty.  One  simple  and  commonly  used  represen¬ 
tation  of  this  effect  is  the  modified  form  of  the  formula  for  the  burning  rate 
given  by  (3-5).  Under  unsteady  conditions ,  visual  observations,  particularly 
of  metallized  propellants,  suggest  that  the  behavior  is  likely  to  much  more 
complicated.  Consequently,  the  simple  response  function  based  on  un¬ 
steady  heat  conduction  in  the  solid  is  probably  not  an  accurate  representation 
of  the  behavior.  Additional  time  lags  associated  with  processes  at  the 
interface  between  the  solid  and  the  gas,  and  possibly  in  the  gas  phase  itself, 
are  probably  im.portant. 
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Unfortunately,  no  quantlative  Information  about  those  processes  is 
available.  For  the  purposes  here,  then,  a  very  simplified  approach  is  taken. 

It  is  assumed  that  the  transient  behavior  is  due  solely  to  the  unsteady  thermal 
conduction  in  the  solid  phase.  While  this  approximation  may  not  adequately 
mo'  el  the  total  response  of  real  propellants,  there  Is  little  /eason  to  doubt 
that  the  thermal  wave  must  be  present  in  any  case.  An  assessment  of  the 
importance  of  the  thermal  wave  response  relative  to  the  complete  velocity 
coupled  response  of  a  propellant  must  await  further  study. 

This  does  not  mean  that  the  response  to  velocity  coupling  wll’  be  either 
qualitatively  or  quantitatively  similar  to  that  for  pressure  coupling.  Ihe  reason 
for  this  is  that  velocity  coupling  is  intrinsically  non-linear,  a  point  which  has 
been  emphasized  in  Refs.  27-31.  This  non-linearity  arises  because  the  burn¬ 
ing  is  sensitive  to  the  magnitude  but  not  the  direction  of  the  velocity  parallel 
to  the  surface.  It  Is  therefore  independent  of  details  of  the  physical  processes 
and  is  always  present.  This  is  essentially  the  only  feature  peculiar  to  velocity 
coupling  which  has  been  discussed  quantitatively  in  the  literature.  In  the  pre¬ 
sent  work,  the  intent  is  to  obtain  numerical  results  for  comparison  with  results 
obtained  for  pressure  coupling.  The  latter  is  ot  course  treated  here  as  a  purely 
linear  phenomenon. 


In  view  of  the  above  remarks,  and  with  appeal  to  the  discussions  of 
Refs.  27-31,  the  response  of  the  mass  flux  to  transient  disturbances  of  velocity 
parallel  to  the  surface  is  written,  by  analogy  with  the  result  (4-34)  for  pressure 
coupling,  as 


X 

X 


X+V  -  (1+A)+AB, 


•[ej(jul-u^)  -  e2(  u  -u^)  ]•  (4-71) 


Here  u^  is  the  threshold  velocity.  If  has  been  introduced  in  previous 
work  as  a  result  of  experimental  observations  on  steady  burning  which 
have  shown  that  the  effect  of  erosion  on  the  burning  rate  appears  in  some 
cases  to  be  absent  unless  the  parallel  velocity  u  is  larger  than  some 
value  u^.  Ihe  quantities  and  62  are  defined  as  follows: 


(4-72) 


fO  |ul<u 
®1“^1  lul>u^j 
rO  |u|<u^  . 

~  ^1  iuj>u^  J 

The  subscript  ()v  on  n  and  B  is  used  to  Indicate  that  the  values  of  these 
parameters  may  differ  from  those  for  pressure  coupling.  In  the  limit  (U-»o  -  i.e. 
a  step  change  of  u  is  imposed  -  Equation  (4-71)  becomes 

^  =n^6(|ul  -u^)  (4-73) 

m 

The  change  of  course,  is  all  in  juj;  u^  is  held  constant.  If  the  threshold  velocity 
is  introduced  in  (3-5),  the  linearized  form  gives 

G* 

6r  =  6m  =  n  6p  +  k _  6(lu|-u^)  (4-74) 

r  m  p  l+C*(|uj-u^) 


Comparison  of  (4-73)  and  (4-74)  shows  that  the  correct  limit  is  obtained  only  If 


the  parameter  n^  is  given  by 


"v  = 


1  +  C*  (juj  -  u^ 


(4-75) 


It  obvious  that  for  velocity  coupling  alone,  the  formula  (4-68)  can  be  used 
to  compute  the  transient  mass  flux,  but  with  n  replaced  by  n^,  B  by  B^  and 
by  Cj  dul  -  u^)  -  ^2  Oul  -  u^).  ^ 

When  pressure  and  velocity  coupling  are  accounted  for  simultaneously,  and 
they  are  assumed  to  be  additive,  the  transform  of  the  mass  flux  is  given  by 


X  +  Y  -  (1+A)  +  AB  p 


n  AB 

_ V  V _ 

X  +  ^  -  (1+A)  +  AB 

K  V 


/ £j(lul-up  -e2(u-up] 


The  total  transient  burning  response  can  then  be  calculated  as 

—  '  —  'p  '  —  'v 

m  m  m 


(4-76) 


(4-77) 


vvhere  (-^— )  ,  the  pressure  coupled  response,  is  calculated  from  equation  (4-68), 
and  (^^the  velocity  coupled  response  is  similarly  calculated  by  making  the 
propePsubstitutions  (indicated  previously)  in  equation  (4-68). 
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5. 


BOUNDARY  CONDITIONS 


5 . 1  Head  End  -  x  =  0 

The  boundary  conditions  at  x  =  0  assume  the  head  end  is  rigid.  The 
boundary  condition  for  the  gas  phase  is  then, 

X  =  0  u  =  0  (5-1) 

The  “correct"  boundary  condition  for  the  particles  at  the  head  end 
would  have  to  allow  for  particles  colliding  with  the  wall.  The  collision 
could  be  assumed  to  be  completely  elastic,  inelastic,  or  anywhere  in  - 
between.  However,  if  particles  were  allowed  to  bounce  off  the  wall, 
there  would  clearly  exist  two  classes  of  particles  at  locations  near  the 
end  wall;  those  that  were  appro^'ching  the  wall,  and  those  that  had  already 
been  reflected.  The  existence  of  particles  having  different  velocities 
at  the  same  location  is  not  compatible  with  the  assumptions  underlying 
the  two-phase  flow  analysis,  as  outlined  in  Section  3.1,  and,  hence, 
cannot  be  accommodated  in  the  present  model.  With  this  in  mind  the 
particle  boundary  condition  has  been  simply  taken  to  be, 

X  =  0  u  =  0  (5-2) 

This  approximation  should  not  seriously  affect  the  overal’  solution. 

5.2  End  of  the  Grain  -  x  =  1 

The  boundary  condition  at  the  end  of  the  grain  (x  =  1,  in  the  non- 
dimensional  coordinate  system)  is  approximated  as  follows.  In  rocket  motors 
with  short  nozzles  the  convective  derivatives,  u  ,  are  much  larger  than 
the  time  derivatives,  ,  in  the  flow  field  between  the  end  of  the  grain  and 
the  throat.  A  quasi-steady  approximation  for  the  flow  in  this  region  should, 
therefore,  be  a  reasonably  good  approximation. 

While  it  is  possible  to  obtain  an  "exact"  one-dimensional  two-phase 
nozzle  solution,  ^^^\he  '  fractional -lag"  approximation^^^ ' has  been  shown 
to  yield  very  good  results,  and  is  much  simpler.  The  fractional  lag 


analysis  is  based  on  the  usual  two  phase  flow  assumptions  (see  Section 
3.1),  with  one  important  difference.  It  is  assumed  that  the  particle 
velocity  and  temperature  are  constant  fractions  of  the  gas  velocity  and 
temperature,  respectively. 

u*  T*,-T*  1-T 

^  P  s  _ —  Y 

u*  T£  -  T*  1-T  '  ' 

Q  F  g  g 

It  can  then  be  shown  that 

L  =  [1  +  3  Pr  C(l-K)/Kr^  (5-4) 

where  Pr  is  the  gas  Prandtl  number  and  C  is  the  ratio  of  particle  to 
gas  specific  heat. 

The  fractional  lag  assumption  allows  the  one-dimensional  two-phase 
equations  of  motion  to  be  reduced  to  the  one-dimensional  isentropic 
equations  for  a  perfect  gas  with  altered  properties.  If  m  =  puA  and 
m  =  mass  flow  rates  of  the  gas  and  particles,  respectively, 

at  the  end  of  the  grain,  then  the  fractional  lag  analysir  yields  the  following 
for  the  equivalent  isentropic  exponent,  Y. 

7  =  1  +  (y  -  1)  I  (5-5) 


where , 


B  = 


1  +  ^ 
m 

m 

1  +  ^CL 
m 


m 

E  =  1  +  -^  |K[y(1-K)  +K]  +  (Y-DCLBj- 


(5-6) 


(5-7) 
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and  R  is  the  nondimen sional  throat  radius  of  curvature,  R  =  I'f'  order 

to  find  K  from  equation  (5-9)  7  must  be  known.  However,  equation 
(5-5)  shows  that  K  must  be  known  to  find  y  •  The  establishment  of 
compatible  K  and  y  values  requires  iteration . 

Once  the  values  of  K  and  Y  have  been  established  the  Mach  number 
at  the  end  of  the  grain  is  then  found  by  solving  the  usual  isentropic  equation 
for  Mach  number  as  a  function  of  area  ratio. 


i-A  -,z 

[m  j  = 


JL  C-i. 
m2  Ly« 
e 


J}  Y -1 


(5-11) 


^The  fractional  lag  analysis  yields  the  choking  condition,  M  =  1,  i.e. , 
the  Mach  number  of  the  "equivalent”  perfect  gas  is  unity  at  tne  physical 
throat.  The  actual  gas  Mach  number  does  not  reach  unity  until  some 
point  downstream  of  the  threat. 
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The  actual  gas  Mach  number  at  the  end  of  the  grain  is  then 


M  =  —  M 

e  a  ^  e 

e 


(5-12) 


The  boundary  condition  at  x  =  1  is  then,  that  the  gas  Mach  number 
be  equal  to  that  given  by  equation  (5-12).  This  matching  condition  assures 
that  the  flow  will  properly  choke  at  the  throat  (at  least  within  the  confines 
of  the  present  assumptions  of  one-dimensional  quasi-steady  nozzle  flow 
with  constant  fractional  lag).  The  treatment  of  the  nozzle  flow  and  boundary 
condition  at  the  end  of  the  grain,  as  presented,  is  a  simplified  and 
approximate  one;  but  one  which  should  give  a  fair  representation  of  the 
actual  flow.  If,  after  further  study,  it  can  be  deduced  that  a  more  accurate 
treatment  of  the  nozzle  is  required,  it  will  not  be  inordinately  difficult  to 
do  so. 
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6 


INITIAL  CONDITIONS— STEADY  STATE  SOLUTION 


In  order  to  initiate  the  method  of  characteristics  instability  solution, 
the  steady  state  equations  of  motion  must  be  solved.  The  initial  condi¬ 
tions  are  then  generated  by  adding  a  perturbation  to  the  steady  state 
solution . 

6 . 1  Steady  State  Equations  in  Conservative  Form 

The  steady  state  two-phase  analysis  assumes  that  the  mass  added  to 
the  chamber  (from  propellant  burning)  is  at  constant  enthalpy  and  has  no 
axial  momentum.  The  rate  of  particle  mass  addition  is  taken  to  be  a  constant 
fraction  of  the  gas  phase  mass  addition  rate.  In  addition,  the  typical  two- 
phase  flow  assumptions  are  made,  i.e. ,  the  particles  are  spheres  of  a 
single  size  with  uniform  internal  temperature,  the  particles  do  not  interact 
with  each  other  and  are  of  negligible  volume,  etc.  With  these  assumptions 
the  steady  state  equations  may  be  obtained  from  the  time  dependent 
equations  (Section  3.2)  by  removing  the  time  derivatives.  The  gas  phase 
equations  become: 


Continuity: 


Momentum: 


Energy: 


State: 


^*(p*u*A*)  =a'*A*  (6-1) 

~*(p*U*2a^  =  “  '^*d^* 

^  '  o*u*  (CpT*+  1/2  u*^)A*  I  =  a^CpT*  +  F^(u*-u*)A*  +  Q*A* 

(6-3) 

P*=  p*RT*  (6-4) 


PilpiqnpPIPPliPIPfliflpffPIpni'ip^^  ' 


Many  limes  It  is  advantageous  to  work  with  the  equations  of  motion 
in  conservative  form.  Equations  (6-1)  and  (6-3)  are  already  in  conservative 
form  and  equation  (6-2)  may  be  written  in  this  form  as 


dx 


*  (P*+  p*u’^2)a* 


(6-5) 


The  particle  aquations  are  already  in  conservative  form  and  are; 


Continuity: 

^.(p;  uJA*)=a,*pA* 

(6-6) 

Momentum: 

|j.(p;  A»)  -  -  f;  A* 

(6-7) 

Energy: 

1:*  r  -  q;  a*-f;u;a* 

(6-8) 

6*2  Nondimens  tonal  Equations 

Using  the  definitions  presented  earlier  in  Section  3.3,  the  steady  state 
equations  can  be  written  in  the  following  nondimensional  form  (Note:  pressure 
has  been  eliminated  from  equation  (6-2)  using  the  equation  of  state  (6-4). 


(o  uA)  =  a'  A 

(6-9) 

4  (OpUpA)  =  »pA 

(6-10) 

Ij'^pAd  +  vu^)'  =vAFp  +  (pT)  ^ 

(6-11) 

It  (Pp  u^A)  =  -  AFp 

(6-12) 

3x.PuAT+  .j2  +  (Y-l)Fp(Up-u) +QplA 

(6-13) 

“p}  “  [“p® Vp'’5pj* 

(6-14) 
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Equations  (6-9  to  14)  represent  six  simultaneous  first  order  differential 
equations  which  can  be  conveniently  integrated,  since  they  are  of  the  form 


where 


(i  =  1  ...  6) 


(6-lS 


I*  DU  A 

•i 


PA{T  +  r  u^) 


0  u^  A 
P  P 


puA  [T  + 


(6-16] 


0)  A 
P 

YAF  +  p  T  ^ 

P  ^  dx 

-  AF 

P 

[U.  +  (Y-l)Fp(Up-u)  +Qp]A 
[V''-l)FpUp-QlA 


(6-17 


Equations  (6-9  to  14)  in  the  form  of  (6-15)  are  integrated  using  an  Adams- 
Moulton  integration  routine.  At  each  step  the  solution  yields  new  values 
for  the  f's  (6-16) ,  and  the  flow  variables  themselves  must  then  be  calcu¬ 
lated  so  the  g's  (6-17)  can  be  reevaluated,  and  the  integration  continued. 
Wr.h  the  f^  known  the  flow  variables  may  be  calculated  as  follows: 


u 


P 


The  steady  state  equations  are  integrated  from  x  =  0  to  x  =  1, 

(the  end  of  the  grain)  where  the  gas  Mach  number  must  be  matched  to  that 
of  the  nozale  solution  to  insure  the  choked  flow  condition  is  satisfied. 

As  discussed  in  Section  5.2,  the  fractional  lag  two  phase  nozzle 

solution  establishes  the  Mach  number  at  the  end  of  the  grain,  M  ,  once 

© 

the  nozzle  radius  of  curvature  chamber  to  throat  area  ratio,  chamber 
temperature  and  particle  to  gas  weight  flow  ratio  have  been  selected. 
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Successive  solutions  with  varying  values  of  chamber  pressure,  Pp,  are 
then  carried  out  until  the  value  of  predicted  by  the  steady  state  chamber 
solution  matches  the  value  set  by  the  fractional  lag  solution.  In  other 
words,  varying  the  chamber  pressure  varies  the  amount  of  mass  addition 
into  the  chamber  through  the  pressure  dependence  of  the  solid  propellant 
burning  rate  until  the  mass  flow  required  to  choke  the  nozzle  is  achieved. 
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7.  NUMERICAL  SOLUTION 

The  longitudinal  combustion  instability  model  developed  in  the  previous 
sections  has  led  to  a  set  of  nonlinear  partial  differential  equations,  which 
are  coupled  to  the  integral  equation  for  the  transient  burning  response.  In 
general,  closed  form  sol  ;tions  to  this  set  of  equations  cannot  be  obtained 
without  resort  to  further  approximation.  Therefore,  numerical  techniques 
have  been  employed  in  order  to  obtain  instability  solutions  without  additional 
simplification  of  the  model. 

The  numerical  solution  of  the  total  problem  is  made  up  of  several 
distinct  elements;  the  characteristics  equations,  transient  burning  rate 
equation,  and  steady  state  equations.  Each  of  these  is  individually 
discussed,  in  turn,  in  the  remainder  of  this  section. 

7, 1  Method  of  Characteristics 

The  one  dimensional  unsteady  equations  of  motion  are  hyperbolic,  and, 
therefore,  are  amenable  to  solution  via  the  method  of  characteristics.  The 
method  of  characteristics  is  itself,  in  the  broad  sense,  a  standard  numerical 
method,  however,  the  details  of  its  application  to  a  particular  problem  can, 
and  do.  vary  over  a  considerable  range.  Each  application,  and  set  of 
equations,  is  subject  to  its  own  little  innuendos,  which  should  be  taken  into 
account  in  the  formulation  of  a  finite  difference  analog  to  the  equations.  The 
modified  Euler  method  is  the  one  most  frequently  employed  in  method  of 
characteristic  solutions.  Prior  experience  with  problems  of  the  present  type, 
however,  indicates  that  it  is  preferable  to  handle  the  streamline  characteristics 
in  an  implicit  manner.  This  allows  larger  step  sizes  to  be  used,  and  higher 
gradients  can  be  tolerated  without  deleterious  effects. 

The  actual  solution  of  the  .oresent  problem  consists  of  the  repeated 
application  of  three  types  of  calculations: 

1.  Field  Points 

2 .  Right  Boundary  Points 

6.  Left  Boundary  Points 

Each  of  these  calculations  is  described,  in  detail,  below. 
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7.1.1  Field  Points 


m 


A  field  point  is  created  by  the  intersection  of  left  running  and  right 
running  characteristic  lines.  By  definition,  it  cannot  lie  on  one  of  the 
boundaries,  since  a-"  '  ’.oundary  one  of  the  characteristic  families  is  absent. 
Figure  7-1  shows  how  the  finite  difference  mesh  for  a  field  point  is  created. 
Although  higher  order  methods,  involving  more  points,  can,  and  have,  been 
used,  the  usual  character! site  mesh,  as  shown,  utilizes  information  at  two 
adjacent  points  to  solve  for  the  flow  variables  at  a  third  point,  more  advanced 
in  time.  Here,  it  is  assumed  that  all  required  information  is  known  at 
points  1  and  2,  and  the  solution  at  point  3  is  to  be  sought. 


t  A 


X 


Figure  7-1.  Field  Point  Calculation 
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In  the  present  scheme,  the  solution  is  begun  by  assigning  the  flow 
variables  at  point  3  values  equal  to  the  average  of  their  values  at  points  1 
and  2.  Using  equation  (3-32)  for  the  directions  of  the  characteristic  lines, 
the  location  of  point  3  yields. 


^3 


c  tj  -  c  t2  + 
c"^  -  c" 


(7-1) 


Xo  =  X, 


+  c'*‘(t. 


-V 


(7-2) 


where  the  average  slopes  of  the  characteristic  lines  are  (the  plus  and 
minus  signs  are  used  to  denote  right  and  left  running  characteristics, 
respectively) , 


c^  =  [(uj+aj^)  +  (ug  +  a2)]/2 

c”  =  [(Uj  -  aj)  +  (ug  -  a^)]/2 


(7-3) 


Average  values  of  u,  p,  T,  cu  ,  a,  etc.,  along  each  of  the  characteristics, 
are  then  defined,  e.g. , 


u^  =  (uj  +\i^/2 

u~  =  (uj  +U2)/2 


(7-4) 


This  enables  the  compatibility  relations,  equations  (3-33)  and  (3-34)  to 
be  written  as 


(t3  -  tj)  (tg  -  tj) 


=  RHS"^ 


(7-5) 


(P„-P,)  (u„  -  u,) 

O"  ^  ^  ^  - ^  DTJC" 

(tg  -  t2)  (tg  -  t2) 


(7-6) 
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where  RHS^  and  RHS"*  are  the  right  hand  sides  of  equations  (3-33)  and 
(3-34),  respectively,  evaluated  in  terms  of  the  corresponding  average 
values  of  the  variables,  e,g. , 


+  ^[(V-l)(u+-u+)  +a+]  +  ^ 


(7-7) 


and 


G+  =  (Y  p‘^a+)“^ 
G  =  (v  P“  a“)  ^ 


(7-8) 


Equations  (7-5)  and  (7-6)  are  easily  solved  to  yield  the  pressure 
and  velocity  at  point  3 

G'^Pj  +  G*'P2  +  Uj  -  U2  +  RHS'^  At+  +  RHS-At" 


G"^  +  G" 


(7-9) 


U3  =u^  -  G‘^(P3-Pj)  +RHS+At'^ 


(7-10) 


Since  the  burning  rate  w.  Is  only  pressure  and/or  velocity  dependent, 
its  value  at  point  3  could  now  be  calculated,  in  principle.  In  practice, 
the  complexity  of  the  transient  burning  rate  analysis  precludes  it,  from 
an  economic  standpoint.  The  burning  rate  analysis  is  used  to  calculate 
00  only  after  the  completion  of  all  the  characteristics  calculations  for 
two  or  more  time  levels.  When  burning  rate  values  are  calculated  at 
other  times,  they  are  done  so  by  simple  first  order  extrapolation  from 
previous  results. 
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The  temperature  at  point  3  is  found  by  utilizing  the  compatibility 
relation  (equation  (3-37)),  along  the  gas  streamline,  dx/dt  =  u,  (shown 
as  line  segment '3-4  in  Figure  7-1).  The  intersection  of  the  gas  streamline 
passing  through  point  3  and  the  line  segment  joining  points  1  and  2,  shown 
as  point  4,  is  not  known,  a  priori,  and  must  be  located.  The  slope  of  the 
streamline  is  taken  to  be  the  average  of  u^  and  u^,  hence,  the  location  of 
point  4  is  found  to  be 

_r  ^^2  ~  ^1^  ^3  r  ^^2  ~^1^  1 

^4  ~  L^3  ^1  ^  ^1  (x2-Xj)  ~  .5(0^ +u^)  J  L  (X2 -x^  ~.5(u2+u^)J 

(7-:i  1) 

(x .  -  X.) 

^4  “  ^3  .5(u3+u^) 

All  of  the  flow  variables  at  point  4  may  now  be  evaluated  by  lineraly  inter¬ 
polating  between  points  1  and  2.  Average  quantities,  denoted  by  a  0,  are 
defined  based  on  the  values  at  points  3  and  4.  As  mentioned  earlier,  the 
solution  of  the  streamline  characteristics  is  carried  out  "implicitly".  Here  the 
term  "implicit"  is  used  to  denote  the  fact  that  the  temperature  derivative  in 
equation  (3-37)  is  not  evaluated  strictly  in  terms  of  known  quantities.  Rather, 
the  temperature  on  the  right  hand  side  is  approximated  as  (T^  (the  current 
unknown) +T4)/2 .  This  has  the  effect  of  couplii^g  the  value  of  the  derivative 
to  the  value  of  the  variable,  itself,  and  helps  t'  pi'  vent  wild  over  or  under¬ 
prediction.  In  the  manner  discussed  above,  the  -  n  for  T^  is  found  to  be 


T  = 
^3 


[  p-if  U)(t3  -t4)]T^  + 


11:0) 


(P3-P^)  +RHS 


P  +  .  5  U)  (tg  -  t^) 


(7-13) 


*0n  the  first  iteration,  u.  is  set  equal  to  (u, +U2)/2;  a  good  approximation 
when  the  velocity  is  small. 
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■r«-?'r5' >?'"■•■*■' 


I 


where  RHS  is  defined  as 


RHS  =  (t3-t/^  5  |f|j,t-  +  Stl  +  J(v-l)u-=]+(V-l)Fp(ap-a)+Qj 

(7-14) 


The  sound  speed  can  then  be  found  as 


*3  V^S 


(7-15) 


and  the  equation  of  state  yields 


0 

3  T3 


(7-16) 


The  particle  velocity  and  temperature  at  point  3  are  found  from  the 

particle  pathline  characteristic  (equation  (3-38))  and  its  dual  compatibility 

relations  (equations  (3-40)  and  (3-41))  in  a  manner  essentially  identical  to 

the  solution  for  T^.  In  this  case,  the  intersection  of  the  particle  pathline 

through  point  3  and  1-2  ,  shown  as  point  5,  is  located.  Values  at  point  5 

are  established  by  interpolation  and  averages  of  the  variables  at  5  and  3 

are  defined  (denoted  by  a  subscript  av) .  The  solutions  for  u  and  T  were 

^3  "3 

also  done  "implicitly." 

The  equations  for  and  t^  may  be  obtained  by  simply  changing  all  the 
4's  to  5's  in  equations  (7-11)  and  (7-12),  and  ere  not  given  here.  The  final 
equations  for  u  and  T  are  easily  shown  to  be 


u  = 


P  u 

Pav  P5 


_ 5av;  ®5  Pav  ^  5) 


P3 


av 


('3 

^av 


^5^ 


(7-17) 
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It  will  be  recalled  that  lacking  a  state  equation  for  the  particles,  the 

equation  for  p  cannot  be  put  into  characteristic  form.  The  equation  for 
P  * 

p  may  be  written  as  (see  (3-42)) 


6  6 
P  P 


6t 


Bx 


(7-18) 


wherein  the  velocity  derivative  is  treated  as  an  ordinary  variable,  and  is 
placed  on  the  right  hand  side.  In  the  form  shown,  equation  (7-18)  is 
analogous  to  the  compatibility  relations,  equations  (3-40)  and  (3-41),  and 
Pp  may  be  found  in  the  manner  previously  outlined  for  velocity  and  tempe- 
rat^ire,  provided  a  suitable  approximation  can  be  found  for  the  velocity 
derivative. 


Figure  7-2 . 


Calculation  of  Bu  /Bx 
P 


*The  area  derivative  term  has  been  left  out  of  the  numerical  analysis  at  this 
point.  If  varying  area  solutions  are  sought  in  the  future  the  solution  for 
p  ,  shown  here,  will  have  to  be  modified. 

P 


ou 

A  first  order  method  for  calculating  has  been  devised,  as  follows: 
Depending  on  whether  t2  >  or  <  tj  (see  Figure  7-2)  a  new  point,  6,  is  located; 
along  2-3  at  tj ,  if  t^  >  t2 ,  or  along  at  t2 ,  if  t2  >  tj .  The  velocity  at 
point  6  is  found  by  linear  interpolation  and  is 


(ti-t2) 


%6  '  ''P2  (VV  \  ^2 


(t2-tj) 

Pe  Pi  Itg  tj)  P3  Pi  2  1 


(7-19) 


The  location  of  point  6,  itself,  is  given  by 


X5=X2+c-(t,-t2)  tj>t2 


X5=Xj+c’"(t2-t2)  t2>t, 


(7-20) 


and  the  velocity  derivative  can  then  be  approximated  as 


>  u  -  u 

^  _  Pg  Pi 

5x  x^-x,  ^1  ^  ^2 

D  1 


-s  U  -  U 

^  =-f2.  P6  ^ 

?*x  ^2~^6  ^  ^ 


The  density,  p  ,  can  now  be  found  from  equation  (7-18), 


(7-21) 


0^  +  fli!  -i  C 
P-;  I  P^„  ^  Pc 


D  = 

P-a 


L  ^av  *-5  ax 


(7-22) 
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7.1.2  Right  Hand  Boundary 


A  right  hand  boundary  point  is  created  when  a  right  running  character¬ 
istic  reaches  the  boundary  at  the  end  of  the  grain,  x  =  1.  Since  any  left 
running  characteristics  passing  through  such  a  point  would  have  to  come  from 
outside  the  domain  of  interest  the  information  carried  by  them  is  not  available. 
Without  the  compatibility  relation  on  left  running  characteristics,  one  less 
equation  is  available;  this  is  compensated  for  by  the  right  hand  boundary 
condition  given  by  equation  (5-12).  This  boundary  condition  fixes  the  Mach 
number  at  the  end  of  grain,  such  that  the  flow  will  choke  at  the  throat  in  a 
quasi-steady  manner. 
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The  finite  difference  mesh  and  method  of  solution  for  a  right  hand 
boundary  point  are  quite  similar  to  that  of  the  field  point. 


I 


I 

_ I _ ► 

y  =  1  X 

Figure  7-3.  Right  Hand  Boundary  Point  Calculation 


Again,  points  1  and  2  are  known,  and  the  solution  at  point  3  is  sought.  In 
this  case,  x-  -  1,  is  known;  and  u,  =M  a,  must  be  satisfied.  M  is 
known  f.om  the  boundary  condition,  but  a.  is  not  known  unti;'  the  energy 
equation  is  solved  to  yield  T^.  To  begin  the  solution  u^  and  all  of  the  other 
quantities  at  point  3  are  se’  equal  to  the  value  of  the  respective  variables 
at  point  2.  The  location  of  point  3  is  fc.ind  to  be 


X3  =  I 


'3  =■■  '1 


(X3-X1) 


(7-25) 
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where  is  defined  in  equation  (7-3).  Preceding  in  a  manner  similar  to 
that  shown,  previously,  for  the  field  point,  Pg  Is  found  from  the  compatibility 
relation  on  the  right  running  characteristic  1-3. 

P3  =  Pi  +  [RhS‘^(t3-ti)  -  (u3-Ui)3/G‘*‘  (7-26) 

4*  4" 

where  RHS  and  G  have  been  given  by  equations  (7-7)  and  (7-8). 

The  calculation  of  T-,  a-,  p-,  u  .  T  ,  and  p  from  the  streamline 

Pq  P3  Po 

and  particle  pathlines  follows,  identicafly,  the  procedure  outlined  for  the 
field  point  solution.  The  only  additional  step  is  the  calculation  of  U3,  from 
the  boundary  condition,  after  T3  and,  hence,  a,  ,  have  been  found. 

7.1.3  Left  Hand  Boundary 

A  left  hand  boundary  point  is  formed  when  a  left  running  characteristic 
intersects  the  wall,  x  =  0.  Like  the  right  hand  boundary  point  calculation, 
the  relation  lost  through  the  nonexistence  of  one  characteristic  family  (here 
the  right  running  family)  is  replaced  by  a  boundary  condition.  In  this  case, 
the  boundary  condition  is  u  =  0  at  x  =  0.  (See  equation  (5-1)).  As  discussed  in 
Section  5.1,  an  additional  boundary  condition  has  also  been  placed  on  the 
particle  velocity,  i.e. ,  Up  =  0  at  x  =  0.  As  "  result  of  these  two  conditions 
the  line  x  =  0  is  both  a  streamline  and  a  particle  pathline.  The  fixed  location 
of  these  lines,  In  this  case,  makes  the  solution  somewhat  simpler,  since 
the  points  denoted  4  and  5  in  the  field  point  and  right  hand  boundary  point 
calculations  need  not  be  established. 

The  left  hand  boundary  point  calculation  is  initiated  by  setting  the 
variables  at  point  3  equal  to  their  respective  values  at  point  1.  The  location 
of  point  3  is  found  to  be* 


*The  program  currently  allows  the  left  hand  boundary  to  have  a  velocity,  so 
that  piston  problems  can  be  analyzed  (without  particles).  This  option  was 
used  only  for  program  checkout  and  will  not  be  considered  here. 
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Figure  7-4.  Left  Hand  Boundary  Point  Calculation 


I 

1 

where  c  is  defined  in  equation  (7-3).  The  velocity  at  point  3,  u^,  equals  ) 

zero,  by  definition.  The  pressure  can  then  be  found  from  the  compatibility  | 

relation  for  left  running  characteristics  (equation  (3-34)),  and  is 

P3  =  P2  +  [RHS"(t3-t2)  -  U2]/G~  (7-28) 

where  RHS  can  be  found  by  analogy  with  equation  (7-7)  and  G  is  specified 
by  equation  (7-8). 

The  temperature  is  found  by  integrating  the  gas  streamline  compatibility 
relation  (3-37)  from  point  1  to  point  3,  and  is  given  by 


7-i2 


^^3  = 


I  -*^<t3-‘i)1  Ti->-^  (P3-Pl>  (t3-ti)[^  I  f  l;.t-  ^  ^  ^Qq] 


p  +  i  ^(^3-^1) 


(7-29) 


where  ( )  quantities  are  basei  on  conditions  averaged  between  points  1  and 


The  particle  velocity,  .'s  been  assumed  to  be  zero,  hence,  in 

a  similar  manner,  using  equation  i3-41),  T  is  found  to  be 

^3 


i:np-iSp(t3-t,)lIp^-.(t3-t,)[i^-Qp/C] 


»p  i'“p<‘3'*l> 


(7-30) 


The  particle  density  is  found  in  the  general  manner  discussed  in  the 

field  point  calculation;  however,  in  this  case,  a  point  6  does  not  have  to 

be  located,  since  3u  /3x  can  be  simply  approximated  as 

P 


au„  u_ 

The  value  of  can  be  computed  from  equation  (7-22)  by  replacing  the 
5's  with  I's.  dke  the  other  calculations .  the  left  hand  boundary  point 
calculation  is  iterated  until  equations  (7-23)  and  (7-24)  are  satisfied. 


7.1.4  Ordering  of  the  Calculations 

As  mentioned  earlier,  the  method  of  characteristics  solution  is 
obtained  by  the  repeated  application  of  the  three  unit  processes  previously 
■’escribed;  the  field  point,  and  right  and  left  hand  boundary  points.  In 
order  to  achieve  a  useful  solution,  however,  a  logical  method  for  determining 
the  order  in  which  these  calculations  are  to  be  performed  is  required. 

Several  possible  ordering  schemes  are  possible,  each  having  certain 
advantages  and  disadvantages.  The  scheme  outlined  below  was  selected 
as  being  the  most  appropriate  for  the  problem;  all  constraints  considered. 
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After  a  steady  state  solution  has  been  achieved,  an  initial  pertur¬ 
bation  is  added  to  it  to  generate  the  initial  line,  at  t  =  0,  for  the  character¬ 
istics  solution.  The  number  of  points,  M,  used  on  the  initial  line 
determines  the  mesh  spacing  for  the  solution.  The  characteristics  solution 
is  initiated  by  a  field  point  calculation  using  points  1  and  2.  Additional 
field  point  calculations  are  successively  carried  out  until  a  total  of  M-1 
have  been  performed,  the  last  one  involving  points  M-1  and  M.  At  the 
completion  of  this  first  stage  of  the  solution  the  x-t  diagram  has  the  form 
shown  in  Figure  7-5. 


Figure  7-5.  x-t  Diagram  After  First  Time  Step 

Another  series  of  field  point  calculations  is  then  performed  using 
points  1',  2'  N'{see  Figure  7-5).  Upon  completion,  a  left  hand  boundary 
point  calculation  is  carried  out  using  points  1  and  1',  and  a  right  hand  boundary 
point  calculation  using  points  N'and  M  is  also  performed.  At  this  point  the 
x-t  diagram  looks  like  the  one  illustrated  in  Figure  7-6. 


Figure  7-6.  x-t  Diagram  After  Two  Time  Steps  (one  complete 
"cycle"  of  the  characteristics  calculation) 
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The  features  to  be  noted  from  this  figure  are  as  follows: 

1 .  Except  for  the  new  point  M  (right  boundary  point) ,  all  of 
the  other  points  (1-(M~1))  have  about  the  same  value  of 
time. 

2.  At  the  end  of  these  two  rounds  of  calculations  the  number 
of  points  remains  the  same,  however,  the  locations  of  the 
points  are  somewhat  displaced  from  their  original  positions. 

These  features  cause  the  characteristic  mesh  network  to  become  skewed. 
The  skewness  of  the  network  generally  becoming  worse  as  time  passes. 

This  skewing  of  the  mesh  is  a  normal  feature  of  characteristics  solutions, 
and,  as  long  as  it  does  not  become  drastic,  is  not  usually  a  cause  for 
concern.  In  the  present  case,  however,  the  skewness  is  undesirable  due 
to  the  coupling  between  the  floiv  and  the  transient  burning  rate  calculation. 

The  burning  rate  analysis  is  one  dimensional  in  that  it  considers  each  x 
location  independent  of  the  others,  and  neglects  longitudinal  energy  transfer. 
At  each  x  location,  however,  the  transient  burning  rate  is  a  function  of  the 
pressure  history  at  that  location,  for  all  time.  In  order  to  carry  out  the 
burning  rate  analysis,  then,  the  pressure  must  be  known  at  specific  x 
locations,  not  at  the  somewhat  random  intervals  that  result  from  the  usual 
characteristics  mesh.  This  situation  could  have  been  relieved  by  essentially 
carrying  along  two  different  meshes;  the  characteristic  mesh,  and  an  ortho¬ 
gonal  one,  obtained  by  interpolation  within  the  characteristic  mesh,  having 
fixed  X  locations,  for  use  in  the  burning  rate  analysis.  Such  an  approach, 
however,  results  in  large  computer  storage  requirements  and  an  unnecessarily 
complex  code . 


t 


I 


The  remedy  for  this  problem  that  was  adopted  does  not  require  two 
separate  meshes  to  be  carried  along.  Rather,  when  the  point  shown  in 
Figure  7-6  is  reached,  the  mesh  is  rectified  by  interpolation;  thereby  creating 
a  set  of  points  all  at  the  same  time,  and  at  the  same  spacing  as  on  the 
initial  line.  This  pattern  is  then  repeated  until  the  computation  reaches  the 
desired  time  level.  The  details  of  the  interpolation  procedure  are  presented 
in  the  following  s.-ction. 


ll 
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Interpolation  -  Rectification  of  the  Characteristics  Mesh 


The  interpolation  of  the  characteristics  results,  in  order  to  obtain  a 
rectilinear  mesh,  could  be  carried  out  to  varying  degrees  of  accuracy.  It 
was  felt  that  the  increased  complexity  and  computation  time  associated  with 
higher  order  interpolation  schemes  was  not  warranted  since,  in  many  practical 
cases,  their  theoretical  accuracy  advantages  are  often  not  realized.  Thus 
linear  interpolation  has  been  used  throughout.* 


The  interpolation  is  performed  in  the  following  manner.  First,  values 

of  the  variables  at  the  original  x  locations  are  obtained  by  interpolating 

along  lines  connecting  the  new  points  1,  2,  ...  M.  Each  of  these  points 

(denoted  by  primes  in  Figure  7-7)  will  then,  in  general,  differ  only  slightly 

in  their  time  coordinate;  except  for  point  M,  at  the  right  hand  boundary, 

which  usually  lags  considerably.  The  "regular"  points,  i.e.,  1,  2,  ...  M-1, 

are  then  searched  to  determine  the  one  having  the  smallest  value  of  time, 

denoted  t^,„. 
min 


1  2  3  M-1  M 


Figure  7-7.  Interpolation  Diagram  For  Rectifying 
the  Characteristics  Mesh 


*The  interpolation  procedure  can  affect  the  accuracy  of  the  results.  Test 
cases  with,  and  without,  interpolation,  using  a  quasi-steady  burning  rate, 
have  been  computed  and  a  comparison  of  the  results  showed  essentially  no 
difference,  at  t\e  end  of  200  time  steps.  These  cases  did  not  have  steep- 
fronted  waves;  if  they  did,  differences  between  the  solutions  would 
probably  have  been  evident. 
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This  is  usually  point  1.  Interpolation  is  then  performed  on  each  of  the  x  = 
constant  lines  to  obtain  a  set  of  points  at  the  uniform  time,  These 

points  are  denoted  by  double  primes  in  Figure  7-7 .  The  original  differences 
in  the  time  coordinate  are  shown  greatly  magnified  in  Figure  7-7  so  that  the 
procedure  may  be  clearly  illustrated.  The  new  right  hand  boundary  point, 

M",  is  established  by  extrapolation.  This  could  be  avoided  by  solving 
for  an  extra  "dummy"  right  hand  boundary  point  (shown  as  M')  using  points 
M  and  M-1,  and  then  interpolating.  This  latter  procedure  is  probably  some¬ 
what  more  accurate;  however,  it  has  not  yet  been  incorporated  into  the 
computer  program  (and  probably  will  not  be  unless  the  extrapolation  procedure 
is  found  to  be  seriously  lacking) . 

When  the  double  interpolation  procedure  is  completed  the  set  of  points 
obtained  are  directly  analogous  to  the  original  set  of  points  on  the  initial 
line,  t  =  0.  The  whole  computation  cycle  can  then  be  repeated  as  many  times 
as  desired.  In  the  nondimensional  coordinate  system  employed  the  points 
on  the  initial  line  are  a  distance  Ax  =  1/(M-1)  apart.  It  turns  out  that  the 
time  increment  fora  computational  cycle,  in  the  illustration,  is  almost 
exactly  equal  to  Ax.  This  result  can  be  used  as  a  rule  of  thumb  to  determine 
how  many  computation  cycles  are  required  to  jach  a  given  value  of  time. 

It  should  also  be  pointed  out  that  since  the  time  step  is  directly  proportional 
to  the  spatial  increment,  the  total  number  of  mesh  points,  and  the  computation 
time,  for  a  given  problem  varies  as  Ax®  .  (If  the  number  of  points  on  the 
initial  line  is  doubled,  the  computation  time  will  approximately  quadruple). 

As  a  consequence  of  this  situation,  efforts  to  optimize  the  accuracy-mesh 
size  tradeoff  are  well  rewarded. 

Incidentally,  the  interpolation  procedure  also  proi  :es  two  desirable 
by-products.  First,  the  randomness  of  the  usual  skewed  characteristics 
mesh  makes  the  interpretation  of  the  calculated  results  difficult.  In  fact, 
in  many  applications  the  characteristics  solution  is  stored  on  tape  and  later 
processed  by  a  separate  interpolation  program.  Second,  the  interpolation 
procedure  allows  the  mesh  to  advance  equally  in  time  across  the  entire 
domain  of  interest.  This  eliminates  the  requirement  for  additional  logic 
to  prevent  the  calculation  of  points  past  the  desired  maximum  time  in  some 
areas  of  the  mesh,  while  the  remaining  portions  of  the  mesh  catch  up. 
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7.3  Transient  Burning  Rate  Solution 

The  interpolation  procedure,  just  described,  provides  the  pressure  and 
velocity  histories  required  for  the  transiet  burning  rate  calculation.  However, 
before  the  burning  rate  expression  can  be  numberically  evaluated  the  deriva¬ 
tives  in  equation  (4-68)  must  be  written  in  difference  form,  A  simple  two 
step  formula  is  used,  which  yields  the  following* 


_d  ^  P(t)-.P(.t-6) 

dT  'p'  6 


—  (— )  = 
dT  ^m^ 


^t)  -  ^{t-6) 
m  m 


(7-32) 


Since  m'/in(T)  is  the  unknown  being  evaluated,  the  term  involving  it  is 
brought  over  to  modify  the  multiplier  on  the  left  hand  side  of  equation  (4-68). 
The  burning  rate  expression  is  also  written  in  terms  of  the  nondimensional 
time  T  (see  equation  (4-59),  the  t  in  this  equation  is  dimensional),  while 
the  remainder  of  the  analysis  uses  the  nondimensional  time  t  (equation  (3-18)). 
The  two  times  are  related  as  follows: 


T  = 


(7-33) 


Since  r  is  a  function  of  axial  distance  the  relative  time  constant,  t  ,  also 
varies  axially.  Defining, 


6*  =  6/t 


^=2  = 


^r2  ^c 


C,=  2nABt 
3  c 

|nAB(t^/TT)- 


(7-34) 


*The  details  of  the  solution  are  worked  out  only  for  the  pressure  coupled  case. 
The  calculation  of  m.'/m  for  velocity  coupling  is  essentially  Identical  if  n  ,B 
and  ej(iu|-u^)  -  62 3^®  used  Instead  of  n,B  and  p'/p.  ^  ^ 
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and  using  equations  (7-32)  and  (7-33),  the  transient  burning  rate  can  be 
expressed  in  the  final  form  used  on  the  numerical  solution. 


-1 

e•'(t)=[l-C^^^•(2-t/)+CJ6■]  {fCj  I  f(?)  d5 


(7-35) 


t-6' 


-  C, 


t.(l-t) 


d§  +  Cj6'^^(t-6’)  +  C^6'^  |’(t-6'; 


■^f^(t)[c/^(2-t^6')  +  Cjf]} 

The  first  two  integrals  in  (7-35)  must  be  integrated  from  0  to  Ihe  current 
time,  t  (actually  t-6'),  each  time  step.  As  time  goes  on  these  integrals  take 
longer  and  longer  to  do,  and  require  increasing  amounts  of  computer  memory 
to  store  all  of  the  required  past  history. 

The  second  two  integrals  in  equation  (7-35)  may  be  written  as 


with  their  dependence  on  the  parameter  t  explicitly  removed.  Hence,  these 
two  integrals  need  not  be  reevaluated  from  ?  ’  each  time;  rather,  a  running 

sum  is  kept,  which  is  added  to  at  each  time  ,'p.  The  latest  value  of  the 
sum  of  the  two  integrals  is  then  multiplied  by  the  pre-exponential  factor 
involving  t. 
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The  integrals  are  currently  evaluated  using  the  trapezoidal  rule. 
Numerical  tests  have  confirmed  that  this  simplest  of  integration  formulas 
yields  accurate  results  for  most  cases  of  interest.  If  cases  are  encountered 
for  which  this  method  proves  inadequate,  it  would  not  be  too  difficult  to 
replace  it  with  a  more  accurate  (but  more  complex)  formula . 

7.3.1  Extending  the  Transient  Burning  Rate  Solution  to  i.ong  Times 

The  current  method  for  computing  the  pressure  and  velocity  coupled 
transient  burning  response  requires  four  quantities  to  be  stored  at  every  axial 
location  each  t.’me  ‘■hat  a  calculation  is  made.  These  quantities  are: 


P 


/  m'  . 
anu  {-r— ) 

m  ^ 


(7-37) 


€j(lul-up  -  C2(u-up  and  i~)^ 


m 

As  time  Increoses,  computer  storage  requirements  will  eventually  exceed 
the  available  core  of  the  computer  being  utilized.  Additionally,  as  time  in¬ 
creases,  the  burning  rate  integrals  (Eq.  7-35)  must  be  evaluated  over  longer 
and  longer  intervals.  This  causes  the  program  execution  time  to  increase  non- 
llnearly,  i.e,  if  the  final  time  of  the  solution  is  doubled  the  execution  time 
increases  by  a  factor  greater  than  2.  A  method  for  at  least  partially  ameliorating 
these  drawbacks  has  been  developed,  however. 


I 
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In  practice  the  effect  of  a  disturbance  at  a  given  time  does  not  last  for¬ 
ever.  Its  effect  on  the  burning  rate  should  decrease  as  time  goes  on,  eventually 
damping  out.  Thus,  one  should  be  able  to  limit  the  amount  of  'tack  history  that 
need  be  considered.  The  time  interval  over  which  the  effect  of  past  distur¬ 
bances  disappears  is  dependent  upon  the  values  of  A  and  B  in  a  rather  complex 
manner.  Based  on  the  analytical  results,  presented  in  Appendix  A,  for  the 
special  cases  of  step  or  experimential  pressure  changes  it  can  be  hypothesized 
that  the  effect  of  a  disturbance  at  a  given  time,  t  ,  should  decay  exponentially 
with  time  at  a  rate  something  like  e^  ~y"-l)(T  'or^hgj-r 


i 


x  =  A(l-B) 

y  = 
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(7-38) 


Figure  7-8  shows  the  locus  of  the  successive  peaks  in  the  transient  burning  rate 
response  to  a  pressure  disturbance  of  the  form 

=  .1  cos(nt)  ts3.48  (ts.274) 

(7-39) 

=  0  t>3.48  (t>.274) 

P 

This  response  was  calculated  using  a  small  computer  program  which  was  devel¬ 
oped  solely  to  integrate  equation  (4-68)  for  a  given  pressure  variation. 

The  result  shown  was  obtained  with  the  burning  rote  parameters  A  and  B 

equal  to  15  and  .7  respectively.  With  these  values  of  A  and  B  x  =  4.5  and  y  = 

5.454,  therefore,  based  on  the  aforementioned  hypothesis  it  would  be  expected 

that  the  response  would  die  out  at  a  rate  like  The  actual  response, 

shown  in  Figure  7-8,  damped,  essentially  exponentially,  at  a  rate  between 
-9  5t  -IOt 

e  *  and  e  ;  in  relatively  good  agreement  with  expectations. 

The  example  just  discussed  serves  to  Illustrate  that  the  response  to  a 
disturbance  which  occurs  at  a  particular  time  eventually  damps  out.  Since  the 
current  transient  burning  rate  analysis  is  linear,  this  feature  can  be  exploited 
to  limit  the  amount  of  computer  storage  required  to  obtain  solutions  over  long 
times.  By  using  the  principle  of  superposition  the  pressure  or  velocity  distur¬ 
bance  history,  at  o  given  axial  location,  can  be  considered  to  be  the  sum  of  a 
series  of  discrete  disturbances,  each  localized  in  time.  The  response  to  a  dis¬ 
turbance  which  exists  only  from  t  =  t  ^  to  t  =  t2  will  damp  out  (to  within  any 
given  tolerance)  in  some  time  interval  after  t2,  denoted  At^.  The  size  of  this 
interval  being  a  function  of  A  and  B,  Thus  knowledge  of  the  disturbance  which 
existed  between  t^  and  t2  and  the  subsequent  response  to  it  becomes  superfluous 
after  a  time  equal  to  t2  +  At^. 

Based  on  this  concept,  a  method  for  obtaining  instability  solutions  out  to 
as  large  a  time  as  desired  can  be  constructed,  which  requires  no  more  storage 
than  would  ordinarly  be  required  to  compute  the  solution  out  to  t  =  2At^.  The 
procedure  is  as  follows: 
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Locus  Of  Trdnsient  Buining  Rato  Response  i'*eaks 


Figure 


7-3  Example  Of  Transient  Response  Damping  Subsequent  To  Ihe  Cutoff 
Of  The  Pressure  Perturbation. 


7-22 


1.  For  times  less  then  the  burning  response  is  computed  in  the 
straightforward  manner  discussed  earlier.  Ihis  Is  shown  pictoriaily 

in  Figure  7-9,  for  the  case  of  a  pressure  coupled  response.  (Actually, 
the  disturbance  will  typically  have  several  cycles  In  ^t^,  however 
only  one  is  illustrated). 

2.  At  times  between  t  =  and  t  =  the  disturbance  history  at  a 

given  location  is  divided  into  two  discrete  parts  (see  Figure  7-10). 

One  Qi  ’turbance  beginning  at  t  =  0  and  ending  at  t  =  At .,  and  the 

d 

other  beginning  at  t  =  At^.  As  shown  in  Figure  7-3,  after  At^  the 
response  due  to  the  first  disturbance  begins  to  damp  while  that  due 
to  the  second  disturbance  grows.  The  total  transient  response  is 
simply  the  sum  of  the  two,  i.e. 


rjnli  rJELl-i  +  rjnli 

^c;.al  ■'l  ^2 

m  mm 


(7-40) 


3.  When  t  oecomes  greater  than  2At^  the  respr  .se  to  the  first  distur¬ 
bance  should  have  diminished  to  the  point  where  it  is  negligible. 
Iharefore,  the  calculation  of  ft  is  ceased;  the  second  disturbance 
becomes  the  first,  and  a  new  disturbance,  now  considered  the  second, 

is  initiated  at  t  =  2Atj.  Hie  situation  it  tnen  exactly  as  it  was  at 

d 

t  =  At^.  This  process  is  repeated  after  each  At^^  time  interval  until  the 
desired  time  is  reached. 


When  velocity  and  pressure  coupling  are  treated  simultaneously  the  afore¬ 
mentioned  procedure  leads  to  a  total  response  given  by  the  sum  of  four  contribu¬ 
tions  (att>At^),  i.e. 


tf  Ua.  =  *  cf  ^2  * 


in' 


m‘ 


m’ 


(7-41) 


where  the  p  and  v  denote  pressure  coupled  and  velocity  coupled  response,  re¬ 
spectively;  while  the  subscripts  1  and  2  pertain  to  the  first  or  second  distur¬ 
bance  of  each  type. 
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7,4  Steady  State  Solution  :  I 

- -  , 

The  method  of  solving  the  steady  state  equations  of  motion  has  already  j  | 

I  ^ 

been  discussed  in  Section  6.  Basically,  the  procedure  is  as  follows.  The  |  'i 

I 

equations  of  motion,  in  conservative  form,  are  written  as  a  set  of  six  j  | 

simultaneous  first  order  differential  equations  (equation  6-15);  which  are  -| 

ft 

then  integrated  using  either  an  Adams-Moulton  or  Runge-Kutta  integration 

fi 

routine.  t 


As  discussed  in  Section  6,  the  solution  is  actually  an  iterative  one; 
wherein  a  value  of  chamber  pressnie,  which  yields  the  proper  Mach  number 
at  the  end  of  the  grain,  must  be  found.  Newton’s  method  is  used  to  provide 
further  "guesses"  of  the  chamber  pressure,  after  the  initial  solution  has  been 
carried  out.  This  technique  has  worked  quite  well.  Typically,  about  two 
iterations  are  all  that  is  required.  Even  when  very  poor  initial  pressure  guesses 
were  specified,  solutions  rarely  took  more  than  six  iterations  to  converge. 

There  was,  however,  one  problem  with  the  steady  state  solution.  The 
boundary  condition  specified  at  x  =  0,  u  =  Up  =  0,  corresponds  to  an  initial 
state  of  dynamic  equilibrium  between  the  particles  and  gas.  Like  the 
chemical  reaction  rate  equations  these  equations  are  of  the  "stiff"  type 
and  are  very  difficult  to  handle,  at  near  equilibrium  conditions.  Standard, 
explicit,  integration  routines  like  Runge-Kutta  and  Adam-Moulton  are 
unstable  under  such  conditions,  unless  extremely  fine  integration  steps 
are  used.  This  problem  can  be  overcome  by  the  use  of  an  implicit  inte¬ 
gration  technique,  as  proven  by  the  results  of  Reference  67,  for  the  chemical 
kinetics  system  of  equations.  In  the  present  case,  however,  the  numerical 
difficulties  at  near-equilibrium  conditions  have  been  overcome  by  using  an 
approximation,  which  serves  quite  well.  The  time,  and  expense,  of  changing 
to  an  implicit  formulation  and  solution  was  thereby  avoided. 

The  approximate  method  for  overcoming  the  numerical  problem  is  based 
on  the  observation  that  the  gas  velocity  increases  in  a  very  nearly  linear 
manner  in  constant  cross-sectional  area  chambers.  This  is  due  to  the  fact 
that  the  burning  rate  does  not  vary  to  any  great  extent  from  the  beginning 
to  the  end  of  the  grain,  for  most  motors. 
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It  will  be  recalled  (see  Section  5)  that  the  fractional  lag  nozzle  solution 
is  based  on  the  assumption  that  velocity  is  proportional  to  axial  distance  in 
the  region  of  the  throat.  The  fact  that  u  is  also  approximately  proportional 
to  X  in  the  chamber  means  that  a  fractional  lag  analysis,  similar  to  that  of 
the  nozzle  region,  should  yield  a  very  good  approximation.  The  fractional 
lag  analysis  in  the  chamber  differs  somewhat  from  that  in  the  nozzle  due  to 
the  pressure  of  mass  addition,  however,  very  similar  results  are  obtainea. 
Starting  from  the  fractional  lag  assumptions 


t; 

f 


n 


I 

I 


u  =K'u 
P 

T  =  1  -  L‘(l-T) 
P 


it  follows  from  the  continuity  equations  that 


PoO 


(7-42) 


(7-43) 


where  p,  is  the  ratio  of  particle  to  gas  burning  rate  (equal  to  p^,  the 
particle  to  gas  weight  flow  ratio).  Without  going  into  details,  it  can  be 
shown  that  when  the  burning  rate  w  is  then  assumed  to  be  invariant  with 
respect  to  axial  distance  the  constants  K’  and  L'  become, 

K'  =|r-l  +  (l+4/D)^] 

L'  =  l/"l+3  PrC(l-K)/K'' 

where  the  constant,  D,  is  based  on  a  Stokes  flow  analysis  and  is 


(7-44) 


^  4  0^  oJup/L 

m  b 


(7-45) 


This  approximate  chamber  fractional  lag  solution  is  used  in  one  of 
two  ways.  If  the  particle  size  is  tc  smell,  or  other  conditions  are  such, 
that  near  equilibrium  (almost  no  lag)  conditions  exist  throughout  the  chamber, 
then  the  approximation  is  used  to  obtain  the  complete  solution,  from  x  =  0  to 
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X  =  1.  On  the  other  hand,  if  conditions  are  such  that  near  equilibrium  exists 
only  near  the  head  end,  the  approximation  is  used  only  until  the  normalized 

particle  velocity  is  greater  than  .01  (u  >  .01). 

P 

The  accuracy  of  this  approximate  technique  has  been  tested  in  a  few 
cases  by  comparing  the  results  to  ones  obtained  by  using  the  exact  equations 
and  very  small  step  sizes.  The  so-called  "exact"  solutions  required  that 
the  initial  conditions,  at  x  =  0,  bo  slightly  perturbed  out  of  equilibrium. 

Even  without  accounting  for  the  slight  differences  in  initial  conditions,  the 
"approximate"  and  "exact"  results  agreed  to  within  a  fraction  of  a  percent, 
for  the  conditions  considered.  The  steady  state  particle  and  gas  velocities 
calculated  for  one  of  the  test  cases  are  shown  in  Figure  7-11. 

For  these  conditions  the  gas  velocity  satisfies,  almost  exactly,  the 
assumption  of  linear  variation  with  distance;  the  primary  assumption  in  the 
fractional  lag  analysis.  The  steady  state  solution  for  the  case  illustrated 
was  actually  found  two  ways;  by  using  the  fractional  lag  approximation 

for  the  first  8%  of  length,  and  using  it  over  the  entire  length.  Tne  calculated 
velocities  were  identical;  illustrating  the  accuracy  that  can  be  achieved  with 
the  approximate  technique  when  the  assumptions  upon  which  it  is  based  are 
fairly  well  realized. 


Chamber  Pressure  -  575  psi 
Chamber  Diameter  -  3.4" 

Chamber  Length  -  23" 

Area  Ratio  -  4.5 

Particle  Radius  -  lOu 

Weight  Flow  Ratio  -  .37 


Figure  7-11.  Two-Phase  Steady  State  Solution  Gas  and 
Particle  Velocities  Versus  Distance 


8.  ANALYSIS  OF  LINEAR  STABILITY 

Since  the  numerical  calculation  is  valid  for  all  amplitudes,  it  is  useful 
as  a  check  to  compute  independently  the  result  for  linear  stability  of  a  normal 
mode  in  the  chamber.  Eventually  one  finds  a  formula  for  the  growth  rate  of 
an  initial  disturbance  having  the  spacial  distribution  for  a  classical  natural 
mode.  The  procedure  followed  here  is  that  presented  in  References  31  and  52. 

The  calculation  breaks  down  broadly  into  two  parts:  (1)  linearization 
of  the  governing  differential  equations  and  construction  of  an  inhomogeneous 
wave  equation;  (2)  solution  for  the  complex  wave  number  for  each  of  the 
stationary  (harmonic)  natural  modes.  All  flow  variables  are  assumed  to  be 
sums  of  steady,  (  )  ,  and  fluctuating,  (  )',  parts.  The  fluctuations  are 
taken  to  vary  harmonically  in  time;  thus,  for  example,  the  perturbation  of 
pressure  is 

P  VpJo 

where  (p'/i^Q  is  the  initial  value  and  k  is  the  complex  wave  number, 

k  =  =  (a'~  ia)  (8-2) 

a 

Therefore,  if  i.he  growth  constant,  a  ,  is  positive  the  waves  will  grov/.  For  use 
in  checking  the  numerical  results,  consider  the  value  of  p'  after  one  period 
of  the  oscillation,  so  t  =  l/f  =  Then  since  a/u*  is  normally  much  lees 

than  unity, 
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Hence,  the  calculated  value  of  the  growth  constant  is 


A  value  for  a  is  ,  therefore,  ee-^ily  found  from  the  numerical  calculation;  the 
purpose  here  is  to  obtain  a  formula  for  the  linear  approximation  to  the 
problem . 

The  analysis  begins  with  the  complete  inviscid  conse.  •<  .  n  equations, 
given  in  Section  3.2  but  repeated  here  in  slightly  altered  from. 


Conservation  of  Mass 


(Gas  Phase) 

it  (p»jA)  =  U)A 

(8-4) 

Conservation  of  Momentum 

'’If 

(8-5) 

Conservation  of  Enerav 

(uA)  (8-6) 

a^  2  • 

u>A+QpA+u^  U)p. 

*  'Vs  ■  V  “p* 

where  F  and  Q  are  given  by  equations  (3-8)  and  (3-14)  respectively. 

P 

and  Le  =  e  -  e  =  C  (T  -T  )  is  the  difference  between  the  energy  of  the 
O  OS  O  V  OS  O  u  v» 

gas  entering  at  the  surface  and  that  at  the  same  position  in  the  chamber, 

and  e  is  total  energy  of  the  particles  at  the  surface, 
pos 


*In  the  numerical  analysis,  it  has  been  assumed  that  the  flow  coming  in  at  the 
surface  enters  always  with  the  steady  flame  temperature.  Hence  aT'  =  -  T'. 

It  is  also  assumed  that  ep^^  =  ep. 


An  equation  for  the  pressure  can  be  formed  by  adding  RT  times  (8-4)  to 
R/C^  times  (8-6): 

^  (pA)  +  YP^  (uA)  +  uA  =  (a^  +RAT^  +  u^)  iwA 

(8-7) 

*“p  V 

where,  as  in  the  numerical  analysis,  e^  e  is  assumed. 

pos  p 

The  linearized  forms  of  (8-5)  and  (8-7  )  are 

p  ^  (Au')+A^5'  =  -  p  A-^(uu*)+A(F'-u*w)  (8-8) 

^(Ap')+yp-^  (Au')  =  -  uAfJ-’  -  vp*  ^(uA)+AP‘u  (8-9) 


where 


fii  =  a^  (1)'  +  (.  RT'+RAT')  uu  A  Q' 

'“'v  P 


(8-10) 


Here  the  speed  of  so'  nd  is  for  the  gas  only,  a®  =  y  p/  P  • 

The  wave  equation  for  the  pressure  fluctuation  p’  is  constructed 
by  differentiating  (8-9)  with  respect  to  time  and  substituting  (8-8)  for 
D  (Au')/st: 


A  ^ 


=  h 


11 


(8-U) 


*See  footnote  to  previous  page. 
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with 


h 


11 


p  — —  (u  u')  +  ~ 
cix®  a® 


IP'  i 

at  ax  —8  St  A 
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(8-12) 


-P 


Jl 

ax 


(Gu') 


d  In  A 
dx 


%r  ^tA(F'-u’m)] 


The  boundary  conditions  on  p'  (usually  at  x  =  0,L)  are  set  by  solving  the 
mo;T.enijin  equation  (8-8)  for  ap'/a^). 


^  =  -  f  (8-13) 

^  X 

f  =  F  +7  (uu')  -  (F’  -  u'i)  (8-14) 


A  formula  for  the  wavenumber  is  found  by  comparing  the  problem 
governed  by  Equations  (8-11)  and  (8-12),  with  the  unperturbed  problem 
(i.e. ,  hii  =  f  =  0).  As  shown  in  Ref.  (52),  the  result  is 

L  _ 

i  PjAdx  +  LP^fAdx]  }  (8-15) 

0  ° 

where 

L 

^  pj  Adx  (8-16) 

0 


and  p  ,  k  are  the  mode  shape  and  wavenumber  for  the  unperturbed  problem. 

jC  jv 

With  hu  aj.d  f  given  by  (8-12)  and  (8-16),  the  formula  for  k^can  eventually 
be  put  in  the  form 
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The  first  group  of  terms  represents  the  interactions  between  combustion 
and  the  wave  motions  at  the  head  end  (x  =  0)  and  the  lateral  boundary  (the 
integral  over  oj') .  The  Influence  of  the  nozzle  is  represented  by  the  term 
evaluated  at  x  =  L.  An  important  stabilizing  effect  is  represented  by  the 
Integral  over  (dp^/dx)®.  This  is  due  to  the  Inelastic  process  of  accelerating 
the  incoming  gases  to  the  acoustic  speed  (u'  -dp^/dx).  The  next  group  of 
terms  represents  the  effects  of  particles  in  the  gas  phase  -  a  stabilizing 
influence.  For  precise  correspondence  with  the  numerical  analysis.  Op'  is 
the  fluctuation  of  heat  transfer  from  the  particles  to  the  gas.  The  last  term 
(^T*)  is  due  to  non-isentropic  interactions  between  the  wave  motions  and  the 
combustion  processes  at  the  boundary. 

Now  for  a  uniform  port,  the  unperturbed  mode  shape  and  frequency  are 


p^=  cos  (kj^x) 


(8-20) 


(8-21) 


and  the  acoustic  velocity  is 


= 

dx 


- - sin  (k  .  x) 


(8-22) 


is  conventional  to  introduce  an  admittance  function  for  the  exhaust 


nozzle,  defined  as 


a  VP 


(8-23) 


Then  with  u'  =  u  =  0atx  =  0, 


Itakj  u'pj+—  j  + 

La  Q 


(8-24) 


where  is  the  Mach  number  at  the  entrance  to  the  exhaust  nozzle . 

The  source  of  mass,  oi>,  is  defined  so  that  a'A  is  the  local  rate  at  which 
mass  is  added  per  unit  length  of  chamber.  If  q  is  the  perimeter  and  mj^  is 
the  mass  flux  of  gases  inward  from  the  burning  surface,  then  u;A  =  m^^q.  Now 
the  fluctuations  of  mass  flux,  mj^',  is  related,  for  pressure  coupling  only, 
to  the  pressure  fluctuation  according  to  (see  Section  4.1.4): 


A.  =  R 


(8-25) 


Hence,  the  contribution  from  the  lateral  boundary  in  (8-19)  can  be  written 


The  last  integral  in  Equation  (8-19)  will  in  general  involve  both  the 
admittance  and  response  functions  since,  as  shown  in  the  discussion  of 
the  response  to  transient  burning. 


41!  =  -^  '  M.  A,  -  (R,  -  1)  (GENERAL)  (8-29) 

T  vP 

If  the  processes  are  isentropic,  the  fluctuation  of  flame  temperature,  T^', 
equals  the  local  fluctuation  of  chamber  temperature,  so  ^T'  =  T^'  -  T'  =  0. 

But  as  noted  earlier  the  fluctuations  of  flame  temperature  have  not  been 
taken  into  account  in  the  numerical  analysis:  T^'  =0.  Hence  AT'  =  -T', 
the  isentropic  value  associated  with  the  waves  in  the  chamber: 


^  =  -  (y-  1)  (ISOTHERMAL)  (8-30) 

T  vp 
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Now  k®  -  k®  «(u/a)^  -  k^  -  2  i  ak^/a  and  the  imaginar/  part  of  the  last 
equation  gives 


ol  -  .^b  (r)  „  H  .  —  1 

—  T--Ay®b  -“-<An 


C  (u) 
m  I 


V 


k  L 


-  1  " 
“  V  2  ,  y 


(8-34) 


The  factor  C^(ai^  t)  /  (1  +  is  exactly  the  result  shown  in  Equation  (8-39) 

of  Reference  (S3);  that  is,  the  attenuation  constant  due  to  particles  is 


A 

2  '■'p 


(8-35) 


The  factor  1/2  arises  from  the  integration  over  the  mode  shape,  while  the 

remaining  factors  constitute  an  approximation  to  the  particle  damping.  Hence, 

it  is  clear  that  more  accurate  estimates  of  the  particle  damping  are 

accomodated  in  (8-34)  simply  by  defining  as  shown,  and  using  whatever 

numbers  are  available  for  a  : 

P 


\  ^ 


(A, 


(r) 


M„) 


1  A  T 

T  Ofr,  ^ 

2  _ E _ 


(8-36) 
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This  is  the  final  result  for  a  uniform  port.  Recall  that  it  Involves  explicitly 
the  assumption  of  Isothermal  pressure  coupling. 


8. 1  Numerical  Example 

Consider  the  pulsed  motors  tested  by  Aerojet;  the  following  values 
are  used: 


D  =  1.99  in. 

P 

L  =  23. 5  in. 

p  =  1400  psia 

r  =  .436  in. /sec. 

a  =  3830  ft. /sec. 

T  =  6110  ®R 

p  =  YP/a®=  .01655  slugs/ft. ^  .00852  gm/cm.® 

A  =  3.11  in. ^ 

S.  =  TT  D  L  =  147  in.® 
b  p 

p  =  1.5  gm/cm .  ®  (density  of  propellant) 
s 

The  Mach  number  of  the  gases  leaving  the  burning  surfaces  is  computed  from 
the  continuity  relation  applied  to  the  interfacial  conversion  of  solid  to  gas: 


Then  assuming  that  the  particle  and  gas  velocities  are  the  same. 


u  = 


D  r 
s 


ori + 


(8-37) 


-  y 
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For  the  Aerojet  propellant  [i  /  p  —  .37  and 

P 


,-  _  (1.5)  (.436)  _  .  KK  rw 

u  7"7mflT5T~7i  tTTTTtT  4.66  ft. /sec, 


(8-38) 


The  Mach  number  of  the  gas  flow  is  therefore 


“b'lsl!  =-»“22 


(8-39) 


By  applying  conserv'ation  of  mass  to  the  average  flow  in  the  chamber, 
one  finds  for  Mach  number  of  the  flow  entering  the  nozzle. 


M  =  -4-  M.  .0577 
e  A  0 


(8-40) 


For  quasi-steady  behavior  of  the  nozzle,  the  admittance  function  is  real, 
and  equal  to 


A  =  M^  =  ,00577 

n  2  e 


(8-41) 


With  these  numbers,  (8-36)  gives 


T^=-0289(R^^^>-  1.17)-.0635  -5^ 

(Lateral  Burning)  (Nozzle)  (Particles) 


(8-42) 


A  wave  will  grow  if  a  >  0,  or,  if  the  real  part  of  the  respons 
sufficiently  large, 

/ 

A 

R  W  >  3  37  + 

111 

f 

Equation  (8-36)  gives  the  symbolic  result. 


R, 


(r) 


.  A, 


>1.17+2 


.  a  t 
1  +  - 
M  ^  aM 


e  function  is 


This  shows  the  obvious  stabilizing  effect  of  increasing  Mach  number 
of  the  mean  flow. 
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NUMERICAL  RESULTS 


In  order  to  check  out  both  the  analysis  and  the  programming, 
several  simple  test  cases  were  solved.  First,  the  left  hand  boundary 
condition  was  modified  to  correspond  to  that  of  a  piston  undergoing 
constant  acceleration.  A  perfect  gas  solution  was  then  obtained  without 
propellant  burning  or  particles.  The  flow  field  generated  by  the  piston 
motion  is  isentropic,  and  the  solution  at  any  point  in  the  flow  can  be 
found  using  the  Rieman  invariants.  The  numerical  solution  accurately 
reproduced  the  known  flow  field,  thereby  verifying  the  characteristics 
analysis  by  itself. 

The  programming  of  the  transient  burning  analysis  was  checked  by 
comparing  the  transient  burning  rates  computed  by  the  instability  program 
to  those  obtained  from  a  simple  check  program  (see  Appendix  A).  When 
the  pressure  field  calculated  by  the  instability  program,  at  a  given  x 
location,  was  input  to  the  check  program,  the  results  of  the  two  programs 
were  identical,  and  served  to  corroborate  the  numerical  accuracy  of  the 
transient  burning  rate  solution.  (In  Appendix  A  the  results  of  the  check 
program  were  verified  by  comparison  with  theoretical  solutions  for  several 
special  cases). 

The  aforementioned  test  cases  were  basically  designed  to  check  the 

accuracy  of  "he  computer  programming,  and  do  not,  themselves,  contribute 

to  the  assessment  of  the  present  instability  model.  During  the  current 

project  there  was  not  sufficient  time  to  systematically  vary  all  of  the 

solution  parameters,  or  to  try  to  optimize  the  inevitable  accuracy-solution 

time  trade-off.  However,  a  limited  number  of  instability  solutions  have 

been  obtained  and  the  effects  of  varying  many  of  the  solution  parameters, 

over,  at  least,  a  limited  range,  have  been  demonstrated.  Many  of  the 

solutions  were  obtained  fora  solid  rocket  configuration  for  which  experi- 

(23) 

mental  data  was  available  ;  however,  no  attempt  was  made  during  the 
present  effort  to  obtain  quantitative  results.  Rather,  the  solutions  were 
obtained,  and  compared  with  data,  so  that  the  qualitative  validity  of  the 
results  and  the  predicted  trends  could  be  ascertained.  The  results  obtained 
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to  date  are  summarized  in  the  rest  of  this  section.  Except  for  Section  9.6, 
all  of  the  calculations  considered  pre.ssure  coupling  only. 

9 . 1  Waves  In  a  Closed  End  Tube 

Another,  of  what  might  be  termed  a  consistency  check  was  made  by 
modifying  the  right  hand  boundary  condition  to  correspond  to  a  solid  wall. 
This  allows  the  problem  of  wave  propagation  in  a  closed  tube  to  be  studied. 
If  an  external  means  of  driving  the  waves  were  also  added  problems  of  the 
type  discussed  in  Section  2.1  could  be  analyzed.  No  attempt  to  do  this 
was  made  at  this  time.  For  waves  with  amplitudes  low  enough  for  the 
propagation  speed  to  approach  that  of  an  infinitesimal  sound  wave,  there 
should  be  no  distortion  of  ^he  wave  form  as  it  travels  up  and  down  the  tube; 
although,  if  a  loss  mechanism  is  present,  the  wave  will  eventually  damp 
out.  Without  a  driving  mechanism  waves  cannot  grow,  however,  if  they  are 
of  finite  amplitude  their  wave  form  changes  during  propagation,  and  should 
become  steeper. 

Figures  9-1  and  9-2  show  the  pressure  and  velocity  waves  over  one 
cycle  (the  time  for  a  wave  to  travel  from  one  end  to  the  other  and  back) 
for  the  following  conditions:  length,  23";  chamber  pressure,  360  psi; 
chamber  temperature,  61 10°F;  particle  diameter,  lOu,  particle  weight 
fraction,  .37;  '  =  1.2.  The  initial  pressure  perturbation  corresponds  to  a 
first  harmonic  standing  wave  with  normalized  amplitude  (small)  of  .01. 

At  the  frequency  corresponding  to  the  above  conditions  (about  1400  hz) 

particle  damping  is  large,  and  the  amplitude  of  the  wave  diminishes  rapidly. 

At  these  low  amplitudes  the  flow  remains  linear,  as  can  be  seen  from  the 

invariance  of  the  wave  form  and  the  fixed  location  of  the  node  point  (x  =  .5). 

The  amount  of  particle  damping  predicted  by  the  current  program  was  semi- 

(54) 

quantitatively  compared  to  the  results  of  Dehority  .  The 
agreement  between  the  results  was  within  the  degree  of  uncertainty  to  which 
the  comparisons  were  made  (approximately  10%).  The  pressure  and  velocity 
waves  at  approximately  one  half,  one,  and  two  cycles,  for  a  closed  tube 
with  no  particles  and  an  initial  first  harmonic  pressure  disturbance  of  0. 1 
amplitude,  are  shown  in  Figures  9-3  and  9-4,  respectively.  This  is  not  a 
very  high  amplitude  case  that  rapidly  demonstrates  nonlinear  behavior, 
however,  at  the  end  of  two  cycles  noticeable  steepening  of  the  pressure 
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wave  can  be  seen,  and  the  presence  of  a  second  harmonic  is  easily  seen 
from  the  form  of  the  velocity  perturbation. 

9.2  Effect  of  Burning  Rate  Parameters— A  and  B 

The  transient  burning  rate  analysis  developed  in  Section  4  contains 
two  parameters,  A  and  B,  which  govern  the  burning  response  of  a  propellant 
to  an  applied  pressure  disturbance.  (The  thermal  diffusivity,  h  ,  can  almost 
be  considered  as  an  additional  parameter,  since  its  value  for  a  given 
propellant  is  still  subject  to  a  sizeable  degree  of  uncertainty). 

Depending  upon  the  values  of  A  and  B,  the  transient  burning  rate 
perturbation  can  lead  or  lag  the  pressure  disturbance,  with  amplified  or 
diminished  response  (i.e. ,  the  amplitude  of  the  burning  rate  perturbation 
can  be  larger  or  smaller  than  the  amplitude  of  the  pressure  wave).  In  a 
linear  analysis,  with  the  pressure  assumed  to  be  harmonic,  equation  (4-34) 
yields  the  relationship  between  the  burning  rate  and  pressure.  Implicit  in 
this  expression  is  the  assumption  that  the  pressure  has  been  harmonic  for 
all  time;  or,  in  more  physical  terms,  that  the  process  has  been  going  long 
enough  for  the  initial  behavior  associated  with  the  origin  of  the  pressure 
pulse  to  be  inconsequential.  In  the  present  analysis,  the  pressure  variation 
need  not  be  harmonic,  and  the  response  at  times  immediately  following  the 
generation  of  a  disturbance  can  also  be  calculated  (Equation  4-68).  It  is 
still  informative,  however,  to  look  a  bit  more  deeply  at  the  linear,  long 
time,  response;  since  many  of  its  characteristics  carry  over,  at  least 
qualitatively,  to  the  more  general  case. 

Figure  9-5  shows  the  real  part  of  the  response  (normalized  by  the 
value  of  the  steady  state  pressure  exponent,  n)  for  A  =  15  and  B  =  .7,  .9, 
!♦/  (n,.,  =  0)  as  a  function  of  a  nondimensional  frequency,  fi,  defined  as 

r  -  2nfH 
fi - p- 

f,  being  the  frequency  of  the  pressure  wave  in  hz.  Figure  9-6  is  similar, 
but  for  A  =  40,  B  =  .8,  .9,  1.0.  In  Section  4.2  it  is  pointed  out  t''at  A  and 


B  must  satisfy  the  inequality 


A  < 


B+I 

(B-l)^ 


and  that  the  peak  response  occurs  in  the  neighborhood  of  =  aVb.  The 
responses  illustrated  in  Figures  9.5  and  9.G  do  peak  near  the  values  given 
by  this  relation.  Furthermore,  the  magnitude  of  the  peak  value  is  seen  to 
grow  larger  as  the  values  of  A  and  B  come  closer  to  violating  the  above 
inequality.  At  low  frequencies,  to  the  left  of  the  peak,  the  burning 
response  leads  the  pressure,  while  at  high  frequencies  it  lags.  (At  long 
times,  when  equation  (4-34)  is  appropriate). 

The  transient  burning  response  to  more  general  nonlinear  waveforms, 
which  are  not  pure  harmonics,  and  have  variable  frequency,  will  of  course 
be  somewhat  different.  However,  when  faced  with  the  task  of  trying  to  find 
values  of  A  and  B  which  yield  calculated  results  that  correlate  well  with 
engine  firing  data,  knowledge  of  the  burning  response  to  constant  pure 
harmonic  pressure  disturbances  can  be  quite  valuable. 

Time  has  limited  the  extent  to  which  the  effects  of  varying  the  burning 
rate  parameters  could  be  studied,  using  the  present  nonlinear  model.  In 
order  to  make  the  calculations  that  were  performed  more  representative  and 
meaningful  the  engine  geometry  and  propellant  parameters  were  selected 
to  a pproxirre  tely  simulate  the  conditions  of  the  experimental  engine  firings 
reported  in  Reference  23,  In  these  tests  the  engine  was  usually  subjected 
to  pulsing,  however,  the  configuration  selected  was  one  which  was 
spontaneously  unstable,  i.e.,  without  being  pulsed. 

The  main  engine  parameters  were  as  follows: 


Throat  diameter 
Port  diameter 
Grain  length 
Chamber  pressure 
Chamber  temperature 
Burning  rate 

Particle  weight  flow  latio 


0.775  in. 
1.99  in. 
23.5  in. 
1400  psia 
61 lOOR 
.32  (P/500) 
.37 
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In  order  to  bring  out  the  burning  rate  effects  more  vividly  particle 
damping  effects  were  kept  to  a  minimum  by  selecting  a  particle  size  of 
2|a^.  Linear  analysis  (see  Section  8.1)was  used  to  provide  a  guide  as  to 
the  size  of  the  response  factor  needed  to  initiate  instability.  Curves 
such  as  illustrated  in  Figures  9-5  and  9-6  were  then  used  to  select 
values  of  A  and  B  thought  to  correspond  to  stable  and  unstable  response. 

The  A  and  B  values  chosen  were: 

A  =  15  B  =  .7 

A  =  11.5  B  =  .64 

In  order  to  begin  the  calculations  a  steady  state  solution  was  obtained, 
which  was  then,  in  turn,  perturbed  by  the  addition  of  a  first  harmonic 
standing  pressure  wave  with  amplitude  equal  to  one-tenth  of  the  mean 
pressure.  Instability  solutions  were  then  carried  out  for  approximately 
three  cycles  .  Figure  9-7  shows  the  calculated  pressure  histories  at  the 
head  end,  x  =  0,  and  end  of  the  grain,  x  =  1,  for  A  =  15  and  B  =  .7. 

Figure  9-8  contains  the  corresponding  results  for  A  =  11.5,  B  =  .64.  These 
results  correspond  to  the  type  of  measurement  one  obtains  with  a  pressure 
transducer,  except  that  the  values  shown  are  perturbations,  i.e. ,  they 
have  the  unperturbed  steady  state  values  subtracted  out. 

As  expected,  the  results  indicate  that  ifA  =  15,B=.7  the  initial 
(  perturbation  damps.  While,  with  A  =  11.5,  B  =  .64,  the  wave  ultimately 

(  amplifies.  In  the  former  case,  the  wave  damps  approximately  10%  per 

J  cycle,  while  the  unstable  case  shows  an  interesting  nonlinear  phenomenon 


+The  measured  (as  well  as  calculated)  frequency  was  approximately  375  hz. 
I  Reference  54  was  then  used  as  a  guide  in  selecting  the  particle  size, 

5  Other  calculations  were  also  made  with  5u  and  9  m  particles  (see 

I  Section  9.4).  At  the  quoted  frequency,  9u  particles  are  approximately 

I  the  size  corresponding  to  maximum  particle  damping. 

I  "t+Time  has  been  nondimensionalized  by  L/a^;  a,  being  the  gas  only  sound 

I  speed.  If  the  wave  speed  in  the  two  phase  mmure  were  exactly  equal  to 

I  a^  then  a  cycle  would  correspond  to  a  nondimens ional  time  of  2.  Since 

^  the  two  phase  wave  speed  was  less  than  a,  a  cycle  corresponded  to 

I  t  «  2.3. 
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After  the  first  cycle  the  wave  is  almost  unchanged;  however,  the  second 
cycle  sees  the  wave  grow  by  about  40%,  and  by  the  end  of  the  third  cycle 
the  perturbation  has  been  amplified  an  additional  50%.  Figure  9-22  contains* 
results  for  the  same  conditions,  except  for  particle  size,  which  show  that  the 
perturbation  may  actually  first  damp  before  ultimately  amplifying.  It  appears 
that  this  behavior  may  be  related  to  the  "necking"  plip-’cmena  which  can  be 
observed  in  experimental  pressure  traces,  followin  .-.'g.  An  explanation 
for  the  calculated  results,  and,  hopefully,  at  least  ><  .ortial  explanation  for 
the  "necking"  phenomena,  can  be  inferred  from  the  nature  of  the  transient 
burning  response. 

The  normalized  values  of  the  transient  burning  rate  perturbation, 

(as  calculated  using  Equation  4-68),  for  the  same  two  sets  of  A  and  B 
values,  are  presented  in  Figures  9-9  and  9-10.  The  usual  response  factor 
calculations  are  not  applicable  at  short  times  after  the  onset  of  a  finite 
size  disturbance.  They  assume  that  a  harmonic  disturbance  has  been 
present  for  all  time.  Realistically,  ho^vever,  the  propellant  cannot  be 
expected  to  respond  instantaneously  to  a  sudden  pressure  change.  The 
change  in  pressure  alter.s  the  heat  transfer  rate  to  the  propellant  and 
initiates  a  thermal  wave  t'lat  travels  down  through  the  propellant.  The 
initial  response  of  the  propellant,  therefore,  c(;ntains  start  up  "transients" 
which,  if  the  frequency  and  amplitude  of  the  wave  do  not  change  rapidly, 
eventually  die  out.  Only  under  these  conditions  of  constant  or  slowly  varying 
wave  form,  and  damped  initial  transients,  can  the  historic  form  of  response 
function  yield  even  approximately  valid  results. 

As  shown  in  Figures  9-9  and  9-10,  as  well  as  in  Appendix  A,  the 
nature  and  duration  of  the  initial  transient  burning  response  is  a  strong 
function  of  the  values  of  A  and  B,  as  well  as  of  the  nature  of  the  disturbance. 

With  A  =  15  and  B  =  .7,  the  transient  burning  rate  initially  varies  at  a  frequency 
close  to  that  of  the  pressure  wave,  but  with  shifting  phase  angle.  The  ratio  of  the 
magnitude  of  the  burning  rate  fluctuations  to  the  magnitude  of  the  pressure 
perturbation,  | (m'/m)/(P'/P )  { ,  ((m'/m)/(P'/P)  is  simply  related  to  the 
response  function  by  Equation  4-35,  and  will  at  times  be  referred  to  as  the 
response  function)  can  be  scon  (using  Figures  9-7  and  9-9)  to  vary  from 


! 
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a  little  less  than  1,  initially,  to  a  value  of  about  4  at  the  end  of  one  cycle. 
It  then  decreases  over  the  next  two  cycles  to  around  2,  which  is  approxi¬ 
mately  equal  to  the  usual  "long  time"  value  of  this  ratio  for  the  given  A 
and  B  values  and  frequency. 

Here,  the  most  vivid  aspect  of  the  initial  transient,  is  the  peak  in 
response  at  the  end  of  the  first  cycle,  at  a  value  about  twice  that  of  the 
"steady  state"  response.  Even  with  this  initial  boost  in  response  the 
pressure  decays  over  the  first  cycle,  although  not  as  rapidly  as  during 
the  succeeding  cycles.  It  is  worth  pointing  out  here,  that  the  fate  of  the 
pressure  disturbance  depends  not  only  on  the  magnitude  of  the  burning  rate, 
but  also  on  its  phase  relationship  to  the  pressure.  Since  neither  the 
magnitude, nor  the  phase,  of  the  response  remains  constant  it  serves  to 
somewhat  complicate  the  interpretation  of  nonlinear  instability  calculations. 

The  second  set  of  results  presented  in  Figures  9-9  and  9-10 were  for 
A  =  11.5,  B  =  .64.  Comparison  of  the  pressure  and  burning  rate  histories 
at  the  head  end  shows  that  the  response  factor  is  initially  about  1.5,  and 
then  grows  to  around  6.25,  7,7$  and  8.0,  at  the  ends  of  the  first,  second 
and  third  cycles  respectively.  The  "long  time"  value  of  the  response  factor 
for  this  case  is  in  the  neighborhood  of  8,  so  rather  than  exhibiting  an 
initial  peaking  behavior,  these  results  show  a  gradual  climb  to  a  "steady 
state"  value  over  several  cycles.  It  is  this  relatively  slow  build  up  in 
response  that  accounts  for  the  fact  that  the  pressure  does  not  increase 
during  the  first  cycle;  only  to  be  amplified  at  an  increasingly  faster  rate 
during  the  succeeding  cycles.  Although  not  proven  yet,  it  can  be  readily 
hypothesized  thct  the  "necking  phenomena,"  referred  to  earlier,  that  is 
observed  in  pressure  traces  from  pulsed  engines,  can  be  attributed  to  a 
similar  slow  build  up  in  burning  response.  It  is  hoped,  that  in  the  future, 
an  investigation  relating  to  this  phenomena  can  be  carried  out,  and  the 
validity  of  the  hypothesis  proven. 

The  previous  discussion  was  based  on  observation  of  the  pressure 
history  at  a  fixed  location.  It  is  also  illuminating  to  look  at  the  results 
as  a  function  of  axial  distance,  for  different  times.  Figure  9-11  shows  the 
pressure  distribution  in  the  chamber  :ust  after  the  disturbance  was 


initiated  and  at  the  approximate  end  of  each  of  the  first  three  cycles. 

Figure  9-8  demonstrated  that  the  wave  amplitude  was  unchanged  after  one 
cycle.  The  new  figure,  9-11,  shows  that  the  waveform  also  remained 
essentially  the  same  during  this  interval.  At  the  later  times,  t  =  4.6  and 
6.9,  the  wave  is  seen  to  have  grown,  while  the  null  point  moves  towards 
the  left  and  the  positive  portion  of  the  wave  becomes  somewhat  larger 
in  magnitude  than  the  negative  part  '  l^x=ll^’  wave  is  still 

predominantly  a  first  harmonic,  but  it  is  not  pure  any  longei .  Nonlinear 
gasdynamic  effects  have  caused  the  wave  to  steepen  a  bit;  which,  in 
other  words,  says  the  wave  contains  some  higher  harmonic  content. 

Since  the  current  analysis  is  nonlinear,  the  wave  is,  and  must  be, 
treated  in  its  entirety;  however,  it  is  still  constructive,  at  times,  to 
consider  the  wave  as  made  up  of  a  basic  frequency  and  harmonics. 

Consider  Figure  9-12,  which  shows  the  velocity  perturbation  in  the  chamber 
(u-u^_q)  at  the  same^ond  of  cyclo,times.  If  the  pressure  disturbance  was 
a  true  standing  wave  the  velocity  wave  would  be  passing  through  the  null 
position  at  the  end  of  a  cycle.  As  the  number  of  cycles  increase  it  can 
be  seen  that  the  velocity  wave  develops  a  small  remainder  at  the  end  of 
a  cycle  (note  the  scale  factor  on  the  velocity  axis),  which  by  the  end  of 
the  third  cycle  is  beginning  to  look  reminiscent  of  a  second  harmonic 
velocity  wave  form. 

Remembering  the  indications  of  at  least  a  small  amount  c:  higher 
harmonic  content  in  both  the  pressure  and  velocity  waves,  it  is  instructive  to 
look  at  the  wave  form  of  the  mass  burning  rate  perturbation  (jl,  - 
shown  in  Figure  9-13.  The  wave  forms  depicted  are  definitely  purer  than 
those  corresponding  to  the  velocity  and  pressure.  Note  the  central  nodal 
point  at  X  =  .5  has  hardly  been  displaced  even  after  three  complete  cycles. 
Obviously,  the  burning  rate,  for  this  case,  does  not  rapidly  respond  to 
small,  higher  frequency,  disturbances.  This  is  the  first  evidence  that  the 
current  transient  burning  rate  model  yields  a  response  which,  essentially, 

+The  so-called  transient  burning  rate  perturbation  shown  in  Figures  9-5 
and  9-10  is  actually  (u,  -  . 


et?^^-'-| 


is  able  to  discriminate  between  harmonics.  Something  which  it  must  do  if 
it  is  to  be  a  viable  model. 

The  next  section  will  deal  in  further  detail  v/ith  the  response  of  the 
propellant  to  different  harmonics. 


9.3  Eflect  of  Initial  Disturbance  Wave  Form 

Since  neither  the  form  of  the  disturbances  which  occur  randomly  i  i 
a  rocket  engine  chamber,  nor  the  exact  conditions  generated  by  an  experi¬ 
mental  pulsing  apparatus  are  really  known,  it  would  be  accommodating  if 
the  instability  results  calculated  using  the  present  model  were  insensitive 
to  the  nature  of  the  initial  disturbance.  As  it  turns  out,  however,  this  is 
not  the  case;  at  le^i?  over  the  first  few  cycles. 

Admittedly,  the  exact  form  of  a  small  random  disturbance  should 
hardly  matter  in  an  engine  which  is  spontaneously  unstable,  since  the 
tendency  of  an  engine-propellant  combination  will  be  to  amplify  only  those 
modes  which  are  ultimately  observed.  If  an  analytical  model  is  to  be 
truly  representative,  it  too  must  possess  this  properly.  In  such  a  case 
one  should  use  some  general  wave  form,  made  up  of  arbitrary  percentages 
of  the  first  harmonic  and  various  overtones,  as  the  initial  disturbance  and 
let  the  model  dictate  which  mode  or  modes  will  be  amplified  and  which  will 
disappear.  If  a  single  pure  harmonic  is  used  (1st,  or  2nd,  or  3rd,  etc.) 
as  the  initial  disturbance  one  runs  the  risk  or.  not  being  able  to  simulate 
the  actual  response,  since  all  the  rest  of  the  harmonics  will  not  be  {>resent 
initially  and  may  not  be  created,  ev'en  through  nonlinear  modification  of  the 
original  wave  form. 

Similar  considerations  also  hold  with  respect  to  predicting  the  results 
of  pulse  tests.  However,  there  are  further  complicating  factors  in  this  case. 
(Neglecting,  for  the  present,  the  fact  that  pulse  guns  usually  introduce 
foreign  gases).  First,  the  pressure  pulse  i.ntroduced  is  accompanied  by 
velocity  and  temperature  disturbances.  Accurate  representations  of  each 
of  these  is,  of  course,  not  feasible.  Second,  the  location  and  magnitude 
of  the  pulse  can  influence  the  measured  results,  even  after  many  cycles, 
and  to  the  extent  of  even  determining  whether  the  pulse  will  eventually  grow 
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or  decay.  Hopefully,  the  present  model,  with  suitable  modification  , 
if  necessary,  will  be  able  to  serve  the  dual  role  of  both  predicting  and 
helping  to  understand  the  basic  mechanism  underlying  e-ich  of  these 
effects. 

Several  instability  solutions  have  been  oblained  in  an  initial  ■■ffort 
to  understand  the  ciiiracteristic  behavior  of  the  present  model,  as  it  relates 
to  t'le  foreyoing  discussion.  The  engine  go  >metry  and  propellant  properties 
were  the  same  as  those  listed  in  the  previous  section,  and  the  same  two 
pairs  of  A  and  B  values  were  used. 

The  first  solution  to  ho  discussed  was  identical  to  the  A  =  11.5, 

B  =  .64  case  presented  in  the  previous  section,  except  the  initial 
disturbance  vvas  taken  to  be  a  10%  second  harmonic,  instead  of  the  10% 
first  harmonic  employed  e^irlier.  The  pressure  and  transient  burning  rate 
perturbation  histories  at  the  head  end  and  end  of  the  grain  are  presented 
in  Figures  9-14  and  9-15  .  Like  the  earlier  cases,  this  solution  was  also 
carried  out  for  3  cycles,  however,  since  the  frequency  is  twice  as  large, 
the  elapsed  time  is  only  half  ;.hat  of  the  first  harmonic  cases. 

In  compaHng  these  results  for  the  second  harmonic  to  the  earlier 
results  (Figure^  S  and  9-10  one  is  immediately  struck  by  the  fact  that  the 
second  harmonic  damps,  while  the  first  was  amplified.  This  difference  can 
be  traced  to  the  behavior  of  the  transient  burning  response.  At  the  higher 
frequency  the  response  factor  (ratio  of  lormalized  burning  rate  perturbation 
to  normalized  pressure  perturbation)  peaks  in  the  neighborhood  of  2, 
before  decreasing  to  order  unity;  while,  it  will  be  recalled,  that  the 
response  factor  for  the  first  harmonic  rose  steadily  to  a  value  around  8. 
Additionally,  the  response  at  the  primary  frequency  was  approximately 
in  phase,  while  the  second  harmonic  response  exhibits  a  significant  phase 
lag.  This  behavior  v/as,  of  c.iurse,  to  be  expected,  in  view  of  the 
characteristics  of  ;he  usual  linear  response  Tnctor;  as  shown  in  Figures  9-1 
and  9-2.  For  the  engine  configuration  under  e /nsideration ,  the  frequency 


4. 

For  instance,  the  addition  of  velocity  coupling  to  the  transient  burning 
\  rate  model. 

i 
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of  the  first  harmonic  is  very  near  the  peak  in  the  response  versus 
frequency  curve.  The  frequency  corresponding  to  the  second  harmonic 
is  then  significantly  displaced  to  the  right  of  the  peak,  and,  hence, 
results  in  a  much  lower  burning  response.  The  peak  also  corresponds 
to  zero  phase  shift,  while  the  amount  of  phase  lag  rnonotonically  increases 
as  the  frequency  becomes  higher.  Under  these  conditions,  each  of  the 
higher  harmonics  should  be  increasingly  damped,  as  their  frequencies 
move  further  and  further  away  from  the  peak  on  the  response  curve. 


The  pressure  and  velocity  waveforms  in  the  chamber,  for  the  second 
harmonic  case,  are  presented  in  Figures  9-16  and  9-17,  respectively. 

Figure  9-16  shows  the  initial  pressure  disturbance,  and  the  waveform 
at  the  end  of  two  cycles  (t  =  2.  39).  The  velocity  peaks  one-quarter  cycle 
after  the  pressure;  therefore,  in  order  to  show  the  velocity  waveforms 
at  close  to  their  peak  values,  they  are  shown  in  Figure  9-17  at  t  =  .3, 
and  t  =  2.7,  Earlier,  in  discussing  the  type  cf  initial  disturbance  to  select 
when  solving  a  spontaneous  instability  problem,  it  was  suggested  that  a 
general  waveform  containing  the  first  and  several  higher  harmonics  be 
employed.  In  order  to  see  how  this  idea  would  work  out  in  practice,  an 
instability  solution  was  obtained  starting  with  an  initial  pulse  consisting 
of  the  sum  of  a  5%  first  harmonic  and  a  5%  second  harmonic.  The  same  A  and 
B  values  (11.5  and  .64)  and  engine  configuration  were  stipulated.  It  will 
be  recalled  that  individually  the  first  harmonic  was  amplified,  while  the 
second  was  damped.  The  results  presented  in  Figure  9-18  should,  therefore, 
come  as  no  surprise.  This  figure  shows  the  initial  pressure  distribution  in 
the  chamber,  and  its  subsequent  development.  The  curves  are  spaced  in 
time  at  approximate  intervals  of  it  =  1.3,  corresponding  to  the  time  for  a 
full  cycle  at  the  frequency  of  the  second  harmonic.  Every  other  curve, 
then,  corresponds  to  a  full  cycle  of  the  first  harmonic.  It  can  be  seen, 
that  as  time  progresses,  the  presence  of  a  second  harmonic  becomes  less 
and  less  obvious,  until,  at  the  last  time  (t  =  6.90),  it  has  all  but 
disappeared.  The  same  effect  can  be  seen  in  the  velocity,  as  demonstrated 
by  Figure  9-19.  Here,  as  in  Figure  9-12,the  velocity  waves  are  shown 
only  slightly  displaced  from  their  null  positions  (the  times  plotted  are  not 
exactly  full  cycles,  and  there  is  also  some  effect  of  nonlinearities)  so 


9-11 


their  magnitude  is  small.  In  fact,  a  comparison  of  the  pressure  and 
velocity  wave  forms  at  t  «  6.9,  with  those  of  the  corresponding  first 
harmonic  only  (initial  pulse)  case  in  Figures  9-11  and  9-12,  indicates  that 
there  is  less  second  harmonic  content  in  the  waves  which  started  with 
50%  second  harmonic.  This  seeming  paradox  is  due  to  the  presence  of 
nonlinear  effects,  as  explained  below. 

Both  the  "first  harmonic  only"  and  the  "mixed"  case  were  initiated 
with  10%  magnitude  pulses.  In  the  former  case  the  pressure  grov/s,  after 
the  first  cycle,  and  by  t  =  6.9  has  been  amplified  about  220%.  In  the 
"mixed"  case  the  first  harmonic  continues  to  behave  in  the  same  manner, 
however,  the  second  harmonic  content  (50%  initially)  is  almost  completely 
damped  out  by  t  =  6.9.  As  a  result,  the  wave,  after  initially  being  damped, 
has  only  amplified  by  about  17%  at  t  =  6.9.  The  nonlinear  effect  of  wave 
steepening,  is,  of  course,  amplitude  dependent.  The  wave  steepening 
process  corresponds  to  the  generation  of  higher  harmonics.  At  the  higher 
amplitude  of  the  "first  harmonic"  case  second  harmonic  content  is  being 
produced  more  rapidly  than  in  the  lower  amplitude  "mixed"  case.  Thus, 
while  the  second  and  higher  harmonics  damp  fairly  rapidly,  they  will  always 
be  present  to  some  degree  in  high  amplitude  waves,  for  they  are  continuously 
being  generated. 

The  pressure  and  burning  rate  histories  at  x  =  0,  and  1,  are  presented 
in  Figures  9- 20 and  9-21 ,  so  they  may  be  compared  to  the  earlier  ones  for 
the  "pure"  first  and  second  harmonic  cases  (Figures  9-8,  9-10,9-14  and 
9-lc).  The  disappearance  of  the  second  harmonic  content  and  lack  of 
burning  response  to  it,  is  also  vividly  portrayed  in  Figures  9-20  and  9-21. 

It  was  gratifying  to  observe  this  ability  of  the  model  to  selectively 
amplify  the  first  harmonic,  since  the  experimental  pressure  traces  recorded 
for  this  engine  configuration  showed  it  to  be  spontaneously  unstable;  with 
very  little  higher  harmonic  content  in  evidence,  even  after  many  cycles. 

The  idea  of  using  an  initial  pulse  containing  all  the  modes  of  interest 
(at  least  the  primary  and  lower  harmonics,  the  higher  frequency  elements 


9-12 


will  always  be  created  by  nonlinear  effects)  has  proven  to  be  an 
interesting  one,  and  is  probably  worth  adopting  as  standard.  The  effort 
to  perform  further  quantitative  comparisons  with  data  is  certainly 
warranted.  This  will  involve,  among  other  things,  continuing  the 
calculations  for  many  more  cycles  so  that  calculated  and  measured 
growth  rates  and  limiting  amplitudes  can  be  compared. 


9.4  Effect  of  Particle  Size 

For  a  given  particle  weight  fraction,  the  amount  of  damping  is  a 

strong  function  of  both  frequency  and  particle  diameter.  Several  linear 

analyses  of  acoustic  wave  particle  damping  have  been  performed , e . g . ,  Ref.  53, 

These  show  that,  at  a  given  frequency,  the  amount  of  damping  peaks  for 

a  certain  diameter  particle,  and  falls  off,  usually  fairly  raj^.dly,  for  larger 

(69) 

and  smaller  particles.  Nonlinear  damping  calculations  '  have  shown 
that  nonlinear  effects  tend  to  shift  the  particle  diameter  corresponding  to 
maximum  damping  without  affecting  the  maximum  amount  of  damping,  itself. 

When  all  of  the  other  gains  and  losses  are  fairly  well  balanced  the 
amount  of  particle  damping  can  be  the  factor  which  determines  the 
stability  of  the  rocket  motor.  Also,  i£  the  amount  of  particle  damping  is 
large,  it  can  improve  the  stability  characteristics  of  an  engine  even  if  the 
other  gains  are  much  greater  thf  the  losses. 

The  present  analysis  represents  a  step  forward  since  nonlinear  particle 
damping  is  included  in  a  total  instability  model  so  its  effect  is  completely 
coupled  to  the  other  combustion  and  flow  phenomena .  Currently,  the  model 
assumes  the  particles  are  all  of  uniform  size.  This  assumption  was  made 
to  enable  the  basic  model  to  be  assessed  in  as  straightforward  and  simple 
a  manner  as  possible.  There  is  no  theoretical  restriction  of  the  type  of 
particle  size  distribution  which  could  be  incorporated  into  the  instability 
analysis,  e.g. ,  an  analytical  particle  distribution  function,  or  an 
approximation  to  a  distribution  using  a  discrete  set  of  particle  sizes.  The 
use  of  even  two  particle  sizes  would  be  a  significant  step  towards 
realistically  matching  the  bimodal  type  particle  size  distributions 
typically  found  in  solid  rocket  motors. 
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The  instability  calculations  presented  previously  were  all  based  on 
a  particle  diameter  of  2^.  A  few  instability  solutions  were  also  performed, 
varying  only  particle  diameter,  in  order  to  assess  the  behavior  of  the 
solution,  with  respect  to  this  parameter. 

Figure  9-22  shows  a  comparison  of  the  pressure  histrries  at  the 
head  end  for  cases  with  particle  diameters  of  2,  5  and  9  microns.  The 
results  of  Reference  54  indicate  that,  at  the  frequency  and  weight  flow 
ratio  of  these  test  cases,  the  maximum  damping  occurs  for  about  9u 
particles.  The  results  shown  do  indicate  that  as  the  particle  size  goes 
from  2  to  9  microns  damping  does  increase.  The  peaked  nature  of  the 
damping  versus  particle  size  curve  is  not  illustrated,  since  larger  than 
9u  particles  were  not  considered  in  this  brief  initial  look  at  particle 
damping.  With  A  and  B  values  of  11.5  and  .64;  linear  analysis  (see 
Section  8.1)  predicts  that  the  engine  would  be  unstable  for  all  three 
particle  sizes  considered.  As  seen  in  Figure  9-22,  the  present  calcu¬ 
lations  corroborate  the  linear  results.  Thus,  this  particular  engine,  could 
not  be  stabilized  by  efforts  to  vary  aluminum  particle  size.  Increasing 
the  weight  fraction  of  aluminum  would,  of  course,  eventually  provide 
enough  damping,  assuming  no  other  deleterious  effects  were  also  induced. 

The  actual  particle  size  distribution  in  the  experimental  engine  that  this 
case  was  based  on,  is  not  known.  The  engine  firings  were,  however, 
spontaneously  unstable;  thus,  there  is  qualitative  agreement  between  the 
current  findings  and  experiment,  regardless  of  the  particle  size. 

Particle  damping  is,  of  course,  duo  to  the  inertia  of  the  solid  particles, 
which  causes  them  to  exert  a  drag  force  on  the  gas  in  their  efforts  to  keep 
up  with  it.  In  a  constant  velocity  flow  the  particles  would  eventually 
attain  the  gas  speed.  When  the  gas  is  accelerating,  however,  the 
particles  can  never  quite  reach  the  gas  velocity,  no  matter  how  small 
the  particles  or  the  rate  of  acceleration.  Figures  9-23  to  9-25  show  the 
particle  and  gas  total  velocities  (not  velocity  perturbations)  about  one- 
quarter  cycle  past  the  end  of  the  first  cycle  (t  »=  2.9),  whe’’  the  gas 
velocity  is  passing  through  its  peak;  for  each  of  the  three  particle  sizes 
considered.  Figure  9-23  shows  that  the  2^  particles  are  small  enough  to 
just  about  keep  up  with  the  gas.  The  difference  in  particle  and  gas 


velocity  being  only  a  fraction  of  a  percent.  As  the  particles  become 
larger  their  inertia  causes  them  to  lag  further  and  further  behind  the  gas. 
Figures  9-24  and  9-25  show  approximately  a  6%  lag  for  the  5^  particles 
and  a  15%  lag  for  the  9u  particles.  Note  that  the  magnitude  of  the  gas 
velocity  (and  pressure,  see  Figure  9-22)  decreases  as  the 
increased  from  2  to  9  microns,  due  to  the  larger  amount  of  energy 
dissipated  through  drag. 

9.5  Effect  of  Initial  Disturbance  AmpliLude  and  Drag  Law 

The  extent  to  which  pulse  amplitude  serves  as  a  determining  factor 
in  the  stability  of  a  solid  propellant  engine  is  not  yet  understood.  Nor  is 
the  role  of  the  initial  pulse  known  with  respect  to  its  influence  upon  the 
limiting  amplitude.  The  foregoing  are,  of  course,  nonlinear  phenomena, 
and,  as  such,  are  not  amenable  to  solution  by  the  linear  techniques  which 
have  been  the  backbone  of  solid  rocket  instability  analysis,  in  the  past. 
The  present  computer  program  represents  a  tool,  which,  hopefully,  will 
allow  these  phenomena  to  be  investigated  and  the  extent  of  their  influence 
to  be  determined. 

All  that  has  been  done  to  the  present  time,  is  to  look  at  the 
differences  in  response  (over  the  first  cycle  only)  to  initial  disturbances  of 
varying  amplitude.  This  preliminary  effort  serves  only  as  a  relative 
indicator  of  the  extent  to  which  nonlinear  effects  are  encountered  at 
various  amplitudes.  Some  of  these  results  are  presented  below. 

Figure  9-26  shows  a  comparison  between  the  head  end  pressure 
histories  for  two  instability  solutions  which  differ  only  in  the  magnitude 
of  the  initial  pulse.  The  pressure  perturbations  have  been  normalized  by 
these  initial  values,  in  each  case,  to  allow  the  differences  to  be  more 
easily  discerned.  One  solution  has  been  discussed  earlier;  the  10% 
magnitude,  A  =  11.5,  B=.64,  first  harmonic  case.  The  second  solution 
employed  an  initial  2%  pulse.  As  illustrated  in  Figure  9-26,  there  is  not 
too  much  difference  betw  en  the  two  solutions.  The  normalized  ampli¬ 
tudes  of  the  two  solutions  are  essentially  the  same  at  the  end  of  the  first 
cycle,  however,  the  history  of  the  larger  amplitude  wave  shows  it  has 
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become  somewhat  steeper  than  the  smaller  amplitude  wave.  The  differences 
between  the  two  solutions  would  continue  to  grow  as  time  increased,  and 
the  waves  amplified  (this  case  was  unstable).  These  results  do  indicate, 
however,  that  waves  in  the  0  to  10%  amplitude  range  are  not  rapidly 
distorted  due  to  the  influence  on  nonlinear  effects . 

A  similar  pressure  plot,  which  covers  a  wider  amplitude  range,  is 
presented  in  Figure  9-27,  for  somewhat  different  conditions^  The 
chamber  pressure  was  360  psi,  port  to  throat  area  ratio  was  4.5,  and  lOu 
particles  were  used.  The  transient  burning  rate  parameters,  A  and  B,  were 
taken  to  be  10  and  .8,  respectively.  All  other  conditions  were  essentially 
the  same  as  those  for  the  solutions  previously  discussed.  These  new 
conditions  correspond  to  significantly  increased  particle  and  nozzle 
damping  and  decreased  burning  response.  (The  frequency  of  the  first 
harmonic  for  this  case  was  about  1400  hz).  It  is  not  surprising,  therefore, 
that  the  waves  damp  quite  rapidly. 

Two  of  the  curves  shown  in  Figure  9-2  7  correspond  to  first  harmonic 
disturbances  with  initial  amplitudes  of  1%  and  50%  of  the  mean  pressure. 
The  third  curve  also  corresponds  to  a  50%  initial  disturbance,  however, 
in  this  case,  a  Stokes  law  drag  coefficient  was  used  to  see  what  the 
difference  would  be.  (The  drag  law  used  in  all  the  other  solutions  is 
given  by  equation  (3-9)).  It  can  be  seen  that  the  50%  initial  amplitude 
waves  are  subject  to  large  and  rapid  steepening,  even  though  they  have 
damped  to  one-third  their  starting  amplitude  during  the  tirst  cycle.  Although 
not  too  evident  in  the  Figure,  the  50%  initial  amplitude  waves  were  damped, 
percentagewise,  somewhat  more  than  the  1%  wave.  To  be  exact,  the  50% 
waves  were  both  damped  about  67%,  while  the  1%  wave  was  damped  by 
only  63%.  Although  the  conditions  for  which  these  solutions  were 
obtained  are  not  the  most  representative  of  solid  rocket  engines,  in 
general,  the  results  are  interesting  because  they  show  the  current 
model  is  capable  of  predicting  an  amplitude  dependence  for  the  damping 
rate.  This  is  an  important  feature,  one  that  is  necessary,  if  a  model  is 
to  be  capable  of  calculating  limiting  amplitudes. 
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The  results  for  the  two  50%  initial  amplitude  cases  are  also  Interesting 
since  they  illustrate  the  effect  changing  the  drag  law.  At  this  large  amplitude, 
particle  Reynolds  numbers  (based  on  relative  velocity  between  particle  and 
gas)  on  the  order  of  300  are  encountered,  far  beyond  the  range  where  Stokes' 
law  is  valid  (Re  <  1) .  At  such  Reynolds  numbers ,  the  nonlinear  drag  law 
given  by  equation  (3-19)  yields  a  drag  coefficient  several  times  larger  than 
Stokes'  law.  Despite  this  fact,  Figure  9-27  shows  that  even  such  a  large 
change  in  drag  coefficient  did  not  measurably  affect  the  amount  particle 
damping  even  when  particle  damping  was,  itself,  quite  large.  (Particle 
damping  accounts  for  about  2/3  of  the  total  damping  for  this  case).  The 
Stokes  law  wave  is  seen  to  have  steepened  up  considerably  more  than  the 
other  one,  however,  the  amplitudes  at  the  end  of  one  cycle  are  almost 
id  sntical. 

The  relative  insensitivity  of  the  computed  results  to  the  form  of  the 
drag  coefficient-Reynolds  number  relationship  can  be  explained  as  follows. 
Other  things  being  equal,  a  larger  drag  coefficient  will  cause  the  relative 
velocity  between  particles  and  gas  to  decrease.  The  drag  force  exerted  on 
the  gas  is  proportional  to  C^aV^,  therefore,  the  combination  of  higher 
and  lower  uV  creates  at  least  a  partial  balance,  leaving  the  total  force 
relatively  unaffected.  That  the  above  actually  occurs  is  graphically  portrayed 
in  Figure  9-28.  The  solid  lines  denote  the  gas  and  particle  velocities 
calculated  using  Stokes'  law,  while  the  dashed  lines  show  the  same  quantities 
at  the  same  time  as  calculated  using  equation  (3-9).  As  mentioned  previously, 
the  nonlinear  drag  coefficient  is  several  times  larger  than  its  Stokes  flow 
equivalent,  but,  as  seen  in  the  Figure,  the  corresponding  relative  velocities 
have  an  inverse  relationship  to  the  magnitude  of  Cj^. 

In  other  cases,  where  the  waves  grow  or  remain  large  for  may  cycles, 
the  form  of  drag  law  may  be  more  important,  however,  based  on  the  present 
result,  it  does  not  appear  that  the  exact  form  of  the  drag  law  is  critical. 
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9.6  Velocity  Coupling 

A  method  for  calculating  the  nonlinear  effects  of  velocity  coupling  on 
the  growth  of  pressure  waves  in  a  combustion  chamber  has  been  developed 
(see  Section  4.3),  and  a  limited  number  of  solutions  have  been  obtained 
in  order  to  observe  the  qualitative  characteristics  of  the  model.  The  scope 
of  the  present  effort  did  not  allow  for  a  systematic  investigation  of  the 
effects  of  velocity  coupling. 

A  series  of  instability  solutions  have  been  obtained  for  an  engine 
configuration  quite  similar  to  that  described  in  Section  9.2  c  Some  of  the 
dimensions  and  the  propellant  properties  were  modified,  however,  so  the 
results  are  not  directly  comparable.  All  of  the  following  solutions  used 
A  =  11.5,  B=B^=  .64,  n=n^=  .3,  and  2|j  particles .  The  four  solutions 
differed  as  follows: 

1.  Velocity  coupling  only,  u^j^  =  0 

2 .  Pressure  coupling  only 

3.  Pressure  and  velocity  coupling,  u^j^  =  0 

4.  Pressure  and  velocity  coupling,  u^j^  =  .04 

Figures  9-29  and  9-30  show  tl.e  waveforms  of  the  pressure  and  burning 
rate  perturbations,  respectively,  for  the  fourth  case,  i.e. ,  pressure  and 
velocity  coupling  with  u^^  =  .04.  The  four  curves  on  each  graph  depict  the 
results  at  t  =  0  and  approximately  at  the  end  of  each  of  the  first  three  cycles. 

It  can  be  seen  that  the  pressure  initially  damps  during  the  first  cycle,  but 
begins  to  be  amplified  during  the  following  cycles,  as  the  burning  rate  response 
builds  up.  At  the  end  of  the  third  cycle  the  pressure  wave  has  recovered  to 
just  about  its  initial  amplitude.  Looking  at  these  pressure  waveforms  in 
Figure  9-29,  it  can  be  seen  that  very  little  harmonic  distortion  has  been 
introduced  into  the  pressure  wave  over  the  first  three  cycles.  This 
contrasts  with  the  burning  rate  waveforms  in  Figure  9-30,  which  show  a 
marked  degree  of  distortion  (compare  Figure  9-30  with  Figure  9-13). 

The  distortion  in  the  burning  rate  waveform  is,  of  course,  due  to  the 
addition  of  a  velocity  coupled  response  which  is  inherently  nonlinear.  The 
fact  that  the  pressure  wave  had  so  little  harmonic  content  is  evidence  that, 
in  this  case,  the  energy  coupling  between  the  velocity  coupled  response  and 


9-18 


the  pressure  wave  was  quite  inefficient.  Further  evidence  of  this  will  be 
presented  in  later  figures,  as  well  as  an  explanation  of  the  reason  for  this 
weak  coupling.  The  pressure  and  pressure  coupled  burning  rate  histories 
for  this  case  are,  qualitatively,  quite  similar  to  those  shown  previously 
in  Figures  9-8  and  9-10.  Figures  9-31  and  9-32  show  the  time  histories 
of  the  velocity  perturbation  function  (Eqs.  4-71  and  4-72)  and  the  transient 
velocity  coupled  burning  rate  at  several  different  axial  locations.  The 
nonlinear  nature  of  the  velocity  driving  function  is  quite  evident  in 
Figure  9-31.  The  burning  response  tends  to  smooth  out  the  distortions  of 
the  driving  function  as  shown  in  Figure  9-32.  There  remains,  however,  a 
bias  in  the  burning  response  that  results  in  the  waveform  having  a  positive 
contribution  greater  than  the  negative  contribution  during  the  balance  of  the 
cycle.  This  bias  generates  a  net  positive  contribution  to  the  mass  and 
energy  in  the  chamber  and  could  be  responsible  for  the  increase  in  mean 
pressure  observed  in  so  many  instability  traces. 

Figure  9-33  shows  the  pressure  histories  at  x  =  0  for  case  1,  2  and  4. 
The  results  for  case  3  are  quite  close  to  those  of  2  and  4  and  have  been  left 
out  for  the  sake  of  clarity.  Here,  again,  it  is  apparent  that  the  energy 
transfer  from  the  velocity  coupled  response  is  quite  inefficient  in  this  case. 
Thus,  the  pressure  coupled  only  and  pressure  plus  velocity  coupling  pressure 
histories  are  almost  identical;  while  with  velocity  coupling  only,  the  pressure 
damps  rather  quickly  since  there  is  not  enough  gain  to  offset  the  significant 
amount  of  damping  in  the  chamber. 

The  extremely  inefficient  transfer  of  energy  from  the  velocity  coupled 
response  is  due  to  the  phase  relationship  between  the  velocity  coupled  burn 
rate  and  the  pressure.  The  results  of  the  linear  analysis  (Eq.  8-19)  show 
that  the  coupling  between  the  combustion  and  the  wave  motion  is  proportional 
to  (in  nondimens ional  form) 

1 

^  p^  (1)' dx  (9-1) 

0 

If  the  waves  in  the  chamber  were  pure  harmonic  standing  waves  the  pressure, 
velocity  and  velocity  coupled  response  would  vary  as  follows,  in  the 


9-19 


absence  of  rectification  or  threshold  effects. 

«  cos  (k^  x) 

sin  (k^  x)  (9-2) 

(■f'  <'^1  =•> 

I 

Thus,  if  the  velocity  coupled  response  were  linear  (i.e. ,  rectification  and 
threshold  effects  net  .  -esent),  then  Equation  (9-1)  shows  that  the  coupling 
between  the  velocity  coupled  response  and  the  pressure  wave  would  be  zero. 
In  the  present  model  the  amount  of  energy  transferred  to  the  pressure  wave 
thus  depends  upon  how  much  rectification  and  threshold  effects  modify  the 
waveform  of  the  velocity  coupled  response.  For  the  conditions  of  the  present 
solutions  the  coupling  remains  quite  weak.  General  con  sions  about  the 
effect  of  velocity  coupling  must  be  deferred  until  solutions  at  a  range  of 
conditions,  and  for  longer  durations,  have  been  obtained. 
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Figure  9-9.  Transient  Burning  Rate  History  At  Hoad  and  Aft  Ends  of  Grain 
A  =  15,  B  =  .7,  Pn  =  .1  cos  (tt  x),  2  u  Particles 


NORMALIZED  AXIAL  DISTANCE 


Velocity  Waveform  In  The  Chamber 
A  =  11.5,  B  =  .64,  Pq  =  .1  cos  { tr  x ) ,  2  fj  Particles 


ondimonsioncil 


Transient  Burning  Rate  History  At  Head  and  Aft  Ends  of  Grain 
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Figure  9-23.  Comparison  Between  Particle  and 
Gas  Velocity,  t  =  2.9 
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Figure  9-24.  Comparison  Between  Particle  and 
Gas  Velocity,  t  =  2,9 
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Pressure  Waveform 


Velocity  Perturbation  History  at  Several  Axial  Locations, 
A  =  11.5,  B  =  .64,  Pn=.l  cos  (-^  x) ,  2  1 1  Particles 
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10 .  SUMMARY ,  CONCLUSIONS  AND  RECOMMENDATIONS 
10.1  Summary  and  Conclusions 

After  surveying  the  existing  experimental  and  analytical  solid  rocket 
instability  results,  a  new  analytical  instability  model  has  been  developed 
and  solved.  Unlike  previous  analytical  solutions,  the  current  longitudinal 
instability  model  is  nonlinear  and  all  of  the  various  phenomena  affecting  the 
flow  and  stability  of  a  motor  can  be  treated  in  a  coupled  manner,  (All  of  the 
possible  flow  and  combustion  processes  are  not  yet  included,  but  the  model 
is  general  enough'  to  allow  for  their  future  incorporation) , 

The  two  primary  elements  of  the  current  instability  analysis  are  a 
method  of  characteristics  solution  of  the  two  phase  flow  in  the  combustion 
chamber  of  the  motor,  and  a  coupled  calculation  of  a  transient  burning 
rate.  The  transient  burning  rate  analysis  presented,  herein,  is  a  unique  and 
interesting  development.  It  is  based  on  an  extension  of  the  most  popular, 
linear,  harmonic  combustion  response  model.  Ibe  current  method  allows  the 
calculation  of  propellant  burning  response  to  a  pressure  disturbance  of  arbitrary 
waveform,  for  all  time,  including  the  period  immediately  following  the  initiation 
of  the  disturbance.  The  analysis  also  includes  a  model  for  velocity  coupled 
response,  Therefore,  for  the  first  time,  the  nonlinear  effects  of  velocity 
coupling  on  the  growth  of  pressure  waves  in  a  combustion  chamber  can  be  com¬ 
puted. 


The  instability  solution,  itself,  begins  with  the  calculation  of  the  steady 
state  two-phase  flow  in  the  motor.  The  flow  in  the  combustion  chamber  is 
calculated  by  numerically  integrating  the  equations  of  motion,  in  conservative 
form.  The  nozzle  flow  and  choked  flow  condition  are  found  using  the  constant 
fractional  lag  approximation.  The  steady  state  conditions  are  then  perturbed 
and  the  subsequent  wave  motion  in  the  motor  is  calculated  numerically,  using 
the  method  of  characteristics.  The  nature  of  the  engine  response  is  dependent 
upon  the  interaction  the  various  gain  and  loss  mechanisms  in  the  engine;  which 
are,  in  turn,  a  function  of  the  propellant  burning  response,  the  size  and 
amount  of  particulate  matter  present,  the  magnitude  and  shape  of  the  initial 
disturbance  and  the  geometrical  configuration  of  the  motor. 
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In  order  to  proceed  with  the  development  and  assessment  of  the  current 
approach  in  the  most  rational  manner,  a  number  of  simplifying  assumptions 
were  made.  Only  motors  with  cylindrically  perforated  grain  were  considered. 
The  gasdynamic  flow  was  assumed  to  be  one-dimensional  and  the  particles 

in  the  gas  stream  were  taken  to  be  of  uniform  size  and  inert;  thereby  ignoring 
the  processes  leading  to  the  formation  of  the  particulate  matter.  The  nozzle 
Jow  was  assumed  to  be  quasi-steady. 


A  series  of  instability  solutions  have  been  calculated,  wherein  some 
of  the  main  parameters  such  as  particle  size,  burning  rate  constants,  initial 
disturbance  waveform  and  magnitude  and  type  of  response  coupling  have  been 
varied,  in  an  attempt  to  qualitatively  assess  the  behavior  and  validity  of  the 
present  model.  The  results  obtained  are  quite  encouraging.  From  all  appear¬ 
ances,  the  qualitative  behavior  of  the  model  is  quite  realistic  and  comparisons 
with  one  set  of  data  tend  to  corroborate  its  efficacy. 


When  this  effort  was  initiated  it  was  hoped  that  it  would  lead  to  a 
greater  understanding  of  longitudinal  combustion  instability,  and,  thereby 
to  the  realization  of  an  improved  predictive  capability.  It  appears  that  these 
goals  will  be  realized,  since  the  limited  number  of  instability  solutions 
obtained  to  date  have  already  provided  some  new  insights,  or  confirmed  what 
previously  could  only  be  hypothesized. 

For  instance,  when  test  motors  are  pulsed,  to  initiate  an  instability, 
the  initial  disturbance  is  often  observed  to  partially  damp,  before  ultimately 
being  amplified.  Linear  analysis  cannot  account  for  this  behavior,  however, 
as  "necking"  can,  in  all  likelihood,  be  attributed  to  the  existence  of  a  time 
lag  between  the  occurrence  of  a  disturbance  and  the  time  at  which  the  pro¬ 
pellant  burning  rate  can  fully  respond  to  it.  During  this  initial  period  the 
surface  combustion  cannot  supply  energy  to  the  wave  at  a  rate  fast  enough 
to  overcome  the  effect  of  the  always  present  damping  mechanisms.  Given 
time,  however,  the  propellant  response  builds  up,  in  m.any  cases,  to  the 
point  where  a  net  amount  of  energy  is  supplied  to  the  wave,  and  it  amplifies. 
The  duration,  or  even  existence,  of  such  a  response  lag  depends  upon  the 
values  of  the  parameters  appearing  in  the  transient  burning  rate  model.  By 
varying  these  constants  the  initial  response  may  be  made  to  build  up  slowly, 
rapidly,  or  even  overshoot,  the  response  level  it  will  attain  at  later  times. 
Attempts  to  correlate  the  calculated  and  measured  wave  amplitudes  immediately 
following  pulsing  should,  therefore,  provide  some  input  towards  the  empirical 
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determination  of  the  burning  rate  parameters.  The  selection  of  values  for 
A  and  B,  as  these  parameters  are  denoted,  for  quantitative  use,  must  also 
be  based  on  the  ability  of  the  calculated  results  to  match  measured  growth 
rates  and  limiting  amplitudes. 

In  many  instances,  experimental  pressure  traces,  under  unstable 
conditions,  are  observed  to  have  relatively  little  harmonic  content,  even  at 
relatively  high  amplitudes.  This  is  in  contrast  to  the  results  obtained  for 
acoustic  waves  driven  in  closed  resonant  tubes,  where  significant  wave 
steepening,  and  even  weak  shocks,  are  observed  even  at  relatively  low 
amplitudes.  It  has  been  felt  that  this  difference  in  behavior  must  be  attri¬ 
butable  to  some  characteristic  of  the  surface  combustion  response  which 
allows  the  higher  harmonics  to  be  attenuated  while  simultaneously  amplifying 
the  primary  acoustic  mode.  Such  a  cause  and  effect  relationship  has,  however, 
never  been  previously  demonstrated  analytically.  It  was,  therefore,  quite 
gratifying  to  observe  such  an  effect  in  the  present  instability  solutions .  The 
results  presented  earlier  conclusively  demonstrate  that  the  current  transient 
burning  rate  model  is  capable  of  producing  a  response  which  discriminates 
between  the  various  harmonics  contained  in  a  general  pressure  wave; 
amplifying  one,  while  attenuating  several  or  all  of  the  remaining  harmonics. 
These  results,  therefore,  appear  to  confirm  the  hypothesis  that  the  low 
harmonic  content  exhibited  by  many  instability  traces  is  du**  to  the  attenuated 
response  of  the  propellant  to  the  higher  frequency  elements  of  a  disturbance. 

To  date,  instability  solutions  have  not  been  continued  out  to  a  large 
number  of  cycles,  therefore,  the  ability  of  the  current  model  to  predict  a 
limiting  amplitude  has  not  been  demonstrated.  The  present  results  do,  how¬ 
ever,  exhibit  a  nonlinear  feature  which  while  not  guaranteeing  that  a  limit¬ 
ing  amplitude  will  be  reached,  must  be  present  if  such  limiting  is  to  occur, 

Ihel  feature  is  the  amplitude  dependence  of  both  of  the  ga'n  and  damping 
mechanisms. 

It  has  also  been  concluded  from  the  results  that  calculated  particle 
damping  effects  do  not  appear  to  be  sensitive  to  the  exact  form  of  the  drag 
coefficient  relationship.  While  based  on  only  limited  results,  this 
circumstance,  if  corroborated  by  further  study,  ij  a  welcome  one,  since 
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there  is  still  some  disagreement  about  the  proper  variation  ot  with 
Reynolds  number. 

10.2  Recommenda tions 

Based  on  results  achieved  to  date,  the  method  of  analyzing  longi¬ 
tudinal  combustion  instability  developed,  herein,  certainly  merits  further 
investigation;  including  quantitative  comparison  with  experimental  data. 

In  or'^er  to  achieve  such  quantitative  agreement  it  may  be  necessary  to 
elimiuc,  some  of  the  simplifying  assumptions  made  during  this  initial 
effort.  It  may  also  be  desirable  to  relax  some  of  the  other  restrictions, 
currently  imposed,  in  order  to  widen  the  applicability  of  the  method.  Some 
of  the  possible  refinements  or  extensions  of  the  present  work  are  discussed 
below. 

In  order  to  fully  assess  the  current  instability  model,  or  any  further 
modifications  to  it,  solutions  must  be  carried  out  for  many  more  cycles  than 
have  been  computed  to  date.  Only  in  this  way  can  it  be  determined  if  the 
model  is  capable  of  predicting  realistic  growth  rates  and  limiting  amplitudes. 
When  such  solutions  are  sought,  any  improvements  that  could  be  made  in  com 
putational  efficiency  would  be  quite  significant  from  an  economic  standpoint. 
In  this  regard,  it  appears  that  replacement  of  the  method  of  characteristics 
solution  of  the  fluid  dynamic  equations  of  motion  with  a  straight  finite  dif¬ 
ference  solution  of  the  Lax-Wendroff  type  (see  Ref.  68,  for  example)  would 
be  in  order. 

When  the  method  of  characteristics  was  selected  as  the  numerical 
technique  to  be  used  in  solving  the  flow  equations,  the  nature  of  the 
transient  burning  rate  analysis  was,  as  yet,  unknown.  As  it  turned  out, 
the  need  for  pressure  histories  at  fixed  axial  locations  was  satisfied  by 
interpolation,  and  rectification  of  the  characteristics  mesh  at  every  other 
time  step.  These  steps  could  be  elimtinated  if  the  computations  were  per¬ 
formed  in  a  rectilinear  mesh  to  begin  with.  The  difference  solution  of  the 
equations  of  motion  is  also  somewhat  simpler  and  more  straightforward 
with  the  finite  difference  approach.  It  is  estimated  that  a  change  to  a  Lax- 
Wendroff  scheme  of  the  type  discussed  in  Reference  68  could  reduce  the 
computational  time  by  30%  or  more,  from  the  current  level  of  about  2  min¬ 
utes  per  cycle  (  on  a  CDC  6600  computer). 


The  development  of  an  alternative,  or  improved  method  for  calculating 
transient  burning  response  should  also  be  considered,  since  the  relatively  large 
amount  of  computer  storage  currently  required  is  not  a  desirable  feature.  Con¬ 
sideration  should  also  be  given  to  the  development  of  a  nonlinear  transient  re¬ 
sponse  model  for  pressure  coupling.  Further  investigation  of  the  velocity  coupling 
phenomena  should  also  be  considered,  as  the  existing  model  may  not  prove  to  be 
adequate.  A  transient  response  model  based  on  the  nonlinear  solution  of  the  heat 
conduction  equation  may  be  able  to  provide  a  framework  within  which  all  of  the 
aforementioned  items  can  be  accomplished. 

It  may  also  be  possible  to  reduce  the  number  of  axial  locations  at  which 
the  transient  burning  response  need  be  calculated.  Currently,  linear  interpolation 
is  used  to  calculate  burning  rate  at  locations  between  those  at  which  it  is  directly 
computed.  Replacing  the  linear  interpolation  with  a  higher  order  spline  interpolation 
procedure  should  allow  equivalent  accuracy  to  be  achieved  with  larger  spacing  be¬ 
tween  the  locations  at  which  transient  burning  response  is  calculated. 

Other  refinements,  or  extensions,  that  would  be  worthwhile,  and  not 
too  difficult  to  incorporate  into  the  solution  are  the  ability  to  handle  more 
complex  motor  geometries,  and  the  adoption  of  a  more  realistic  particle  size 
distribution.  The  additional  geometries  that  could  still  be  considered  on  a 
one-dimensional  basis  include  more  general  grain  perforations,  and  motor 
cases  having  gaps  in  the  grain,  end  grain,  or  a  grain  that  does  not  extend  the 
full  length  of  the  chamber.  It  would  be  possible  to  incorporate  a  particle  size 
distribution  function  into  the  analysis,  however,  it  is  recommended  that  the 
next  step  in  this  direction  should  be  the  extension  to  two  discrete  particle  sizes 
groups.  Such  an  extension  wouid  be  relatively  simple,  yet  it  would  be  a  sig¬ 
nificant  step  towards  realistically  -.atching  the  bimodal  type  particle  size  dis¬ 
tributions  typically  found  in  sc  li  1  rocket  motors. 

The  development  of  a  model  for  metal  particle  formation  and  burning  is 
not  recommended  at  this  time.  The  whole  mechanism  of  particle  formation 
is  still  subject  to  wide  uncertainties,  and  the  relative  significance  of 
particle  burning  as  a  determinant  of  motor  stability  has  yet  to  be  established. 

If  these  processes  can  be  realistically  modeled  and  are  shown  to  be 
important,  at  some  later  time,  there  is  no  intrinsic  reason  why  they  cannot 
be  incorporated  into  the  instability  model. 


10-5 


‘f^m 


11.  REFERENCES 


6. 


12, 


13, 


Saenger,  R.  A.,  and  Hudson,  G.  E. ,  "Periodic  Shock  Waves 
in  Resonating  Gas  Columns,"  J.  Acoust.  Soc.  of  America, 

V.  32,  No.  8  (Aug.  1960). 


Betchov,  R. ,  "Nonlinear  Oscillations  of  a  Column  of  Gas, " 
Physics  of  Fluids,  V.  1,  No.  3  (May-June,  1958)  pp.  205'-212. 


Chester,  W. ,  "Resonant  Oscillations  in  Closed  Tubes," 
J.  Fluid  Mechanics,  V.  18,  Part  1  (1964)  pp.  44-64. 


Coppens,  A.  B.,  and  Sanders,  J.  V.,  "Finite-Amplitude 
Standing  Waves  in  Rigid-Walled  Tubes,  "  J.  Acoust.  Soc.  of 
Am.  V.  43,  No.  3  (Mar.  1968)  pp.  516-529. 


Weiss,  N.  O.,  "The  Development  of  a  Shock  From  Standing 
Waves  of  Finite  Amplitude  in  an  Isentropic  Fluid,  "  Proc.  Comb. 
Phil.  Soc.  V.  60  (1964)  pp.  129-135. 


Fowl-^r,  J.  R.  and  Rosenthal,  J.  S.,  "Missile  Vibration  Environ¬ 
ment  for  Solid  Propellant  Oscillatory  Burning, "  AIAA  Paper 
71-756  (June  1971). 


7. 


Bergman,  G.  H.  and  lessen,  E.  C.,  "Evaluation  of  Conventional 
Rocket  Motor  Instrumentation  for  Analysis  of  Oscillatory 
Combustion,"  AIAA  Paper  71-755  (June  1971). 


8. 


Browning,  S.  C.,  Kraskin,  M.,  and  Thacker,  J.  H.,  "Application 
of  Combustion  Instability  Technology  to  Solid-Propellant  Motor 
Problems,"  AIAA  Paper  71-758. 


Dickinson,  L,  A.,  Brownlee,  W.  G.,  and  Jackson,  F.,  "CARDE 
Investigations  of  Finite  Wave  Axial  Combustion  Instability," 
Canadian  Armament  Research  and  Development  Establishment, 
Valcartier,  Quebec,  Canada  CARDE  T.  N.  1959/62  (1962). 


10. 


Brownlee,  W.  G. ,  "Nonlinear  Axial  Combustion  Instability  in 
Solid  Propellant  Motors,  "  AIAA  J.  V.  2,  No.  2  (Feb.  1964)  pp.  205- 
284. 


11. 


Brownlee,  W.  G.,  and  Roberts,  A.  K. ,  "Investigations  of 
Nonlinear  Axial  Combustion  Instability  in  Solid  Propellant  Rocket 
Motors,"  CARDE  T.R.  474/64  (June  1964). 


Brownlee,  W.  G. ,  and  Roberts,  A.  K. ,  "Investigations  of 
Nonlinear  Axial  Combustion  Instability  in  Solid  Propellant  Rocket 
Motors,"  AIAA  Paper  63-474  (Oct.  1963). 


Brownlee,  W.  G. ,  and  Kimbell,  G.  H. ,  "Shock  Propagation  in 
Solid-Propellant  Rocket  Combustors, "  AIAA  J.  V.  4.  No.  6 
(June  1966)  pp.  1332-34. 


11-1 


14.  Roberts,  A.  K.,  Bergman,  W.  A.,  and  Brownlee,  W.  G., 
"Longitudinal  Combustion  Instability:  Behavior  of  Propellants 
with  High  Solids  Content,"  CARDE  T.N.  1746/67  (Mar.  1967). 

15.  Kimbell,  G.  H,,  and  Brownlee,  W.  G.,  "Flow  Processes  in  a 
Solid  Propellant  Combustor  During  Unstable  Combustion," 

CARDE  T.N.  1783/68  (Mar.  1968). 

16.  Kimbell,  G.  H.,  and  Brownlee,  W.  G.,  "Status  Report  on 
Optical  Studies  of  Combustion  Instability  at  CARDE", 

CARDE  T.N.  1784/68  (Mar.  1968). 

17.  Roberts,  A.  K. ,  Brownlee,  W.  G.,  and  Jackson,  F., 

"Combustion  Instability  and  the  Design  of  Solid  Propellant 
Rocket  Motors,"  Canadian  Aeronautics  and  Space  Journal, 

V.  16,  No.  1  (Jan.  1970)  pp.  21-27. 

18.  Roberts,  A.  K.,  and  Brownlee,  W.  G.,  "Nonlinear  Longitudinal 
Combustion  Instability:  Influence  of  Propellant  Composition," 
AIAA  J.  V.  9,  No.  1  (Jan.  1971)  pp.  140-147. 

19.  Dickinson,  L.  A.,  "Command  Initiation  of  Finite  Wave  Axial 
Combustion  Instability  in  Solid  Propellant  Rocket  Motors," 

ARS  J.  V.  32  (April  1962)  pp.  643-644. 

20.  Dickinson,  L.  A.,  and  Jackson,  F.,  "Combustion  in  Solid 
Propellant  Rocket  Engines,"  5th  AGARD  Colloquium  on  Combustion 
and  Propulsion  (The  Macmillan  Co.,  N.  Y.,  1963)  pp.  531-550. 

21.  Capener,  E.  L.,  Dickinson,  L.  A.,  and  Kier,  R.  J. ,  "Driving 
Processes  of  Finite  Amplitude  Axial  Mode  Instability  in  Sclid- 
Propellant  Rockets,"  AIAA  J.  V.  5,  No.  5. 

22.  Marxman,  G.  A.,  and  Wooldridge ,  C.  E. ,  "Finite-Amplitude 
Axial  Instability  in  Solid-Rocket  Combustion"  Twelfth  Symposium 
(International)  on  Combustion  (The  Combustion  Institute, 
Pittsburgh,  Pa.,  1969)  pp.  115-127. 

23.  Micheli,  P.  L.,  Private  Communication. 

24.  Price,  E.  W.,  "Axial  Mode,  Intermediate  Frequency  Combustion 
Instability  in  Solid  Propellant  Rocket  Motors,"  AIAA  Preprint 
No.  64-146  (Jan.  1964). 

25.  Angclus,  T.  A.,  "Unstable  Burning  Phenomenon  in  Double-Base 
Propellants,"  Progress  in  Astronautics  and  Rocketry,  V.  1, 

Solid  Propellant  Rocket  Research  (Academic  Press,  1960) 

pp.  527-559. 


11-2 


26.  Watermeier,  L.  A.,  Aungst,  W.  P.,  and  Pfaff,  S.  P.,  "An 
Experimental  Study  of  the  Aluminum  Additive  Role  in  Unstable 
Combustion  of  Solid  Rocket  Propellants,"  Ninth  Symposium 
(International)  on  Combustion  (1963)  pp.  316-325. 

27.  Hart.,  R.  W. ,  Bird,  I.  F.,  and  McClure,  F.  T.,  "The  Influence 
of  Erosive  Burning  on  Acoustic  Instability  in  Solid  Propellant 
Rocket  Motors,"  Progress  in  Astronautics  and  Rocketry,  V.  1, 
Solid  Propellant  Rocket  Research  (Academic  Press,  N.  Y. ,  1960) 
pp.  423-499. 

28.  McClure,  F.  T.,  Bird,  J.  R.,  and  Hart,  R.  W. ,  "Erosive 
Mechanism  for  Nonlinear  Instability  in  the  Axial  Mode  of  Solid 
Propellant  Rocket  Motors,"  ARS  J. ,  V.  32,  No.  3  (Mar  1962) 
pp.  374-378. 

29.  Bird,  J.  F.,  Hart,  R.  W.,  and  McClure,  F.  T.,  "Finite  Acoustic 
Oscillations  and  Erosive  Burning  in  Solid  Fuel  Rockets, "  AIAA  J. , 
V.  3,  No.  12  (Dec.  1965)  pp.  2248-2256. 

30.  Price,  E.  W.,  and  Dehority,  G.  L.,  "Velocity  Coupled  Axial 
Mode  Combustion  Instability  in  Solid  Propellant  Rocket  Motors," 
Proceedings  of  the  2nd  ICRPG/AIAA  Solid  Propulsion  Meeting, 
Anaheim,  California  (1967)  pp.  213-227. 

31.  Culick,  F.  E.  C.,  "Stability  of  Longitudinal  Oscillations  With 
Pressure  and  Velocity  Coupling  in  a  Solid  Propellant  Rocket," 
Combustion  Science  and  Technology,  V.  2,  No.  4  (1970), 

pp.  179-201. 

32.  Sirignano,  W.  A. ,  and  Crocco.  L. ,  "A  Shock  Wave  Model  of 
Unstable  Rocket  Combustors,"  AIAA  J.  Vo.  2,  No.  7  (July  1964) 
pp.  1285-1296. 

33.  Mitchell,  C.  E.,  Crocco,  L.,  and  Sirignano,  W.  A.,  "Nonlinear 
Longitudinal  Instability  in  Rocket  Motors  with  Concentrated 
Combustion,"  Combustion  Science  and  Technology ,  V.  1,  No.  1 
(July  1969)  pp.  35-64. 

34.  Crocco,  L.,  and  Mitchell,  C.  E.,  "Nonlinear  Periodic  Oscilla¬ 
tions  in  Rocket  Motors  with  Distributed  Combustion," 
Combustion  Science  and  Technology,  V.  1,  No.  2  (Sept.  1969) 
pp.  147-169. 

35.  Sirignano,  W.  A. ,  "A  Theory  of  Axial-Mode  Shock-Wave 
Oscillations  in  a  Solid-Rocket  Combustor,"  Twelfth  Symposium 
(International)  on  Combustion  (The  Combustion  Institute, 
Pittsburgh,  Pa.,  1969)  pp.  129-137. 


11-3 


36.  Powell,  ^.A.,  and  Zinn,  B.  T.,  "A  Single  Mode  Approximation 
in  the  Solution  of  Nonlinear  Combustion  Instability  Problems , " 
AIAA  preprint  71-208  (Jan.  1971). 

37.  Powell,  E.  A.,  "Nonlinear  Combustion  Instability  in  Liquid 
Propellant  Rocket  Engines,"  Georgia  Institute  of  Technology , 
Report  GITAER  70-6  (1970). 

38.  Culick,  F.  E.  C.,  "Nonlinear  Growth  and  Limiting  Amplitude 

of  Acoustic  Oscillations  in  Combustion  Chambers,"  Combustion 
Science  and  Techn;,logy,  V.  3  (1971)  pp.  1-16. 

39.  Culick,  F.  E.  C.,  "Nonlinear  Behavior  of  Acoustic  Waves  in 
Combustion  Chambers,"  Work  in  progress,  supported  by 
Hercules,  Inc. 

40.  Priem,  R.  J. ,  and  Guentert,  D.  C.,  "Combustion  Instability 
Limits  Determined  by  a  Nonlinear  Theory  and  a  One-Dimensional 
Model,"  NASA  TND-1409  (1962). 

41.  Weiss,  R.  R.,  and  Osborn,  J.  R.,  "Combustion  Characteristics 
of  Bi-Phase  Rockets,"  AFOSR  70-0768  TR  (1970). 

42.  Hoffman,  R.  J.,  Beltran,  M.  R. ,  Breen,  B.  P. ,  and  Wright,  R. 

D. ,  "Extension  of  the  Priem-Guentert  Annular  Combustion 
Instability  Model  to  a  Bipropeilant  System,"  3rd  ICRPG 
Combustion  Conference,  CPIA  Publ.  No.  138  (1967). 

43.  Beltran,  M.  R.,  et  al,  "Liquid  Rocket  Engine  Combustion 
Instability  Studies,"  Dynamic  Science  Final  Report  AFRPL-TR- 
65-25  (1966). 

44.  Beltran,  M.  R.,  Breen,  B.  P. ,  andGerstein,  M.,  "Liquid 
Rocket  Combustion  Instability  Studies,"  Dvnamic  Science  Final 
Report  No.  SN-68-S1  (1965). 

45.  Priem,  R.  J.,  and  Heidmann,  M.  F.,  "Propellant  Vaporizatioi. 
as  a  Design  Criteria  for  Rocket-Engine  Combustion  Chambers," 
NASA  TR  R-57  (1960). 

46.  Hoffman,  R.  J.,  Wright,  R.  O.,  and  Breen,  B.  P.,  "Combustion 
Instability  Prediction  Using  a  Nonlinear  Bipropellant  Vaporizarion 
Model."  NASACR-920  (1968). 

47.  Campbell,  D.  T.,  and  Chadwick,  W.  D.,  "Combustion  Instability 
Analysis  at  High  Chamber  Pressures,"  AFRPL-TR-67-222 , 
Rocketdyne,  a  Division  of  North  American  Rockwell  Corporation, 
Canoga  Park,  Calif.  (1967). 


11-4 


48.  Berstein,  S.  Z. ,  and  Schechter,  H. ,  "Study  of  High  Frequency 
Nonlinear  Combustion  Instability  with  Baffled  Annular  Liquid 
Propellant  Rocket  Motors,"  Final  Report  MR-7009  (1970). 

49.  Agosta,  V.  D. ,  "Nonlinear  Combustion  Instability:  Longitudinal 
Mode,"  6th  ICRPG  Combustion  Conference,  CPIA  Publ.  No. 

192  (196S). 

50.  Seamans,  T.  F.,  Vanipie,  M.,  and  Agosta,  V.  D.,  "A 
Fundamental  Model  of  Hypergolic  Ignition  in  Space  Ambient 
Engines ,"  AIAA  J. ,  ^,  pp.  1216  (1967). 

51.  Zinn,  B.  T. ,  "A  Theoretical  Study  of  Nonlinear  Traverse 
Combustion  Instability  in  Liquid  Propellant  Rocket  Motors," 
Princeton  University  AMS  iech.  Report  No.  732  (1966). 

52.  Culick,  F.  E.  C.,  "Interactions  Between  the  Flow  Field , 
Combustion,  and  Wave  Motions  in  Rocket  Motors,  "  To  be 
published.  Naval  Weapons  Center,  China  Lake,  Calif.  (1972). 

53.  Dobbins,  R.  A.,  and  Teinkin,  S.,  "Attenuation  and  Dispersion 
of  Sound  by  Particulate- Relaxation  Processes,"  J.  Acoust.  Soc. 
of  Am.,  V.  40,  No.  2  (Feb.  1966)  pp.  317-324. 

54.  Dehority,  G.  L.,  "A  Parametric  Study  of  Particulate  Damping 
Based  on  The  Model  cf  Temkin  ana  Dobbins  Naval  Weapons 
Center,  China  Lake,  California  NWC  TP5002  Sept.  1970. 

55.  Friedly,  J.  C.,  and  Peterson,  E.  E.,  "Influence  of  Combustion 
Parameters  on  Instability  in  Solid  Propellant  Motors.  Part  II, 
Nonlinear  Analysis, "  AI^  J.  V.  4,  No.  11  (Nov.  1966) 

pp.  1932-1937. 

56.  Brown,  R.  H.,  and  Muzzy,  R.  J. ,  "Linear  and  Nonlinear 
Pressure  Coupled  Combustion  Instabilib/  of  Solid  Propellants," 
AIAA  J.  V,  8,  No.  8  (August  1970)  pp.  1492-1500. 

57.  Novozhilov,  B.  V ., "Nonlinear  Oscillations  of  Combustion 
Velocity  of  Powder,"  Soviet  Journal  of  Applied  Mechanics  and 
Technical  Physics,  V.  /,  No.  5,  (May  1966)  pp.  31-41. 

58.  Kriei ,  H.,  et  al,  "Nonstcady  Burning  Phenomenon  of  Solid 
Propellants:  Theory  and  Experiment,"  AIAA  J.  V.  6,  No.  2 
(Feb.  1968)  pp.  278-285. 

59.  Culick,  F.  E.  C. ,  "A  Review  of  Calculations  for  the  Unsteady 
Burning  of  a  Solid  Propellant,"  AIAA  J.,V.  6,  No.  12  (Dec.  1968) 
pp.  2241-2255. 


1.-5 


60.  Culick,  F.  E.  C.,  "Some  Problems  in  the  Unsteady  Burning 
of  Solid  Propellants,"  Naval  Weapons  Center,  China  Lake, 
California,  NWC  TP  4668  (Feb.  1969). 

61.  Price,  E.  W.  "Comments  on  Role  of  Aluminum  in  Suppressing 
Instability  in  Solid  Propellant  Rocket  Motors",  AIA.A  J. , 

Vol.  9,  No.  5,  (May  1971),  pp.  987-990. 

62.  Marble,  F.  E. ,  "Dynamics  of  Dusty  Gases"  Annual  Review 
of  Fluid  Mechanics ,  Vol.  2,  pp.  397-446  Q970). 

63.  Rudinger,  G.,  "Effective  Drag  Coefficients  fcr  Gas-Particle 
Flow  In  Shock  Tubes"  Project  Squid  Technical  Report 
CAL-''?-PU,  March  1969. 

64.  Carlson,  D.  J.  and  Hoglund,  R.  F.,  "Particle  Drag  and 
Heat  Transfer  in  Rocket  Nozzles,  "  AIAA  J. ,  Vol.  2,  Nov. 
1964,  pp.  1980-84. 

65.  Kliegol,  J.  R.  "Gas  Particle  Nozzle  Flows"  Ninth 
International  Symposium  on  Combustion ,  The  Combustion 
Institute,  Pittsburg,  Pa.  pp.  811-826  (1963). 

66.  Kliegel,  I.R.,  "One  Dimensional  Flow  of  a  Gas  Particle 
System,  Paper  No.  60-S,  IAS  28th  Meet'ng,  New  York, 
January  I960. 

67.  Frey,  H.M.,  Nickerson,  G.R.  and  Tyson,  T.J.,  "One 
Dimensional  Kinetic  Nozzle  Analysis  Reference  Computer 
Program,  Dynamic  Science  Report  No.  CS-l-9/’70, 

S'  ..'ember  1970. 

68.  Rubin,  E.L. ,  and  Burstein,  S.Z.  Journal  of  Computational 
Physics,  Vol.  2,  (1967),  pp.  178-96. 

69.  Korman,  H.F.,  "An  Analysis  of  Nonlinear  Particle  Damping" 
unpublished . 

70.  Abramowitz,  M.,  and  Stegun,  I. A.,  (Editors).  Handbook  of 
Mathematical  Functions  with  Formulas  Graphs,  and  Mathe¬ 
matical  Tables .  Dover  Publications ,  Inc. 


11-6 


APPENDIX  A 


THE  TRANSIENT  BURNING  RESPONSE— SOME  SPECIAL  CASES 

The  final  form  of  the  current  transient  burning  rate  model  (Equation 
(4-65),  or  the  approximation  to  it.  Equation  (4-68))  is  not  an  explicit 
formula  for  m'/  in,  and  is  not  very  useful  for  formal  analysis.  For  example, 
it  is  not  possible  to  deduce  from  (4-65)  closed  form  results  for  the  transient 
burning  rate  behavior  at  all  times  t>0.  It  is  relatively  easy,  however,  to 
obtain  short  time  approximations  and  the  long  time  steady-state  behavior,  valid 
after  all  initial  transients  have  died  out.  This  latter  property  may  be  used  to 
advantage  by  comparing  the  long  time  "steady-state"  results  obtained  from 
Equation  (4-65)  to  the  known  response  for  two  special  cases.  Such  a  com¬ 
parison  serves  as,  at  least,  a  partial  consistency  check  on  the  current 
formulation. 

1 .  Response  to  a  Harmonic  Pressure  Oscillation  for  t  w 
For  this  case, 

(A-1) 

P  i 

! 

s 

After  a  long  time,  the  fluctuation  of  mass  flux  also  varies  harmonically,  but  I 

with  phase  shift  i,  and  amplitude  Ir^^!  ; 

! 

Substitution  of  (A-1)  and  (A-2)  into  (4-65)  should  lead  to  the  result  that 
!Rj^|  exp(icp)  is  the  original  response  function  Rj^  given  by  (4-34)  with 


A-1 


r' 


After  the  substitution,  one  haj 


|R,|e''P=  |R^|e‘  [  ^2  e'P?  [  ,9i  ,5] 


(A- 3) 


+  2nAB!  -7 
V 


^  p-^-'  ft  ^ 


S  p"'  «] 


The  integrals  i.iay  be  simply  evaluated,  in  the  limit  as  t  ■*•,  as  follows: 


-tL  '  p~P?  42.  -  _J:_  2 

'^^0  vTvp  - 


V^T 


/;  Vn  ^  i2» 


0 


(A-4) 


e“^'  \  eP"  dv  =  ^  l-e'P”]  .-»  i 


Where  p  =  I  +  i4Q  .  With  the  result  (A-4)  Equation  fA-3)  becomes 


~  J  +  2nAB  1  -pi  +  -^- 


which  can  be  solved  and  rearranged  to  give 


|Rje‘"’=2nAB  - 

P-'"l^p2>i/F*PlP2 


(A-5) 


-  2nAB  7-, 


i  + 


Vp 


(  /p-'C^)(y'p-C'2) 


as  required. 
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2 .  Response  to  a  Step  Change  of  Pressure  for  t  -» eo 
Now  set 


4^  =  A  (t  >  0) 

P 


(A- 6) 


After  a  long  time,  m'/fn  should  approach  a  constant  n  ,  equal  to  nA  .  As 
above,  the  contributions  for  t  small,  from  m'/nT  different  from  p,  are 
supposed  negligible.  Hence,  (4-65)  becomes 


(A-7) 


By  a  simple  transformation, 

T 


0 


and  (A-7)  becomes,  for  large  t 


^  ^  ®  ^  rtco  1 


a  =  p  |o^ +0^  -  0^0^}  +  4nABA 


(A- 8) 


(A-9) 


From  the  definitions  (4-62),  it  follows  that 


“l*  "z  -  ”l'’2  =  ‘  -  "AB 


and 


a  =  p  (1  -4AB)  +  4nABA 


(A- 10) 
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In  addition  to  checking  that  the  numerical  solution  for  the  transient 
burning  rate  aoproached  the  proper  long  time  limits,  indicated  in  Sections 
A.  1  and  A. 2  above,  for  harmonic  and  step  pressure  disturbances,  a  further 
accuracy  check  has  been  made.  The  numerical  results  have  been  compared 
with  exact  results  deduced  from  Equation  (4-61),  for  the  response  to  step 
and  exponential  changes  of  pressure.  All  of  these  comparisons  corroborated 
the  numerical  calculations.  Exact  solutions  for  step  and  exponential 
pressure  changes  may  be  found  as  follows. 


An  exponential  pressure  change  which  approaches  the  finite  value  ti. 
for  long  times  is  represented  by 


£1 

P 


(A-11) 


The  Laplace  transform  in  the  variable  p  is 


00 

\  4e  '  li  (1-e^  ')  e  ^  ’  d' 

0 

(A- 12) 

4a  _  4A 
P-1  p-o^ 

where 

•-^  =  l-e  (A-13) 


Equation  (4-61)  becomes 
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2nAB 


1  + 


’^Vp 


‘•P 


i  1  1  .  D  i'  , 

dp 
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Inversion  is  accomplished  by  using  standard  tables,  after  the  integrand 
has  been  expanded  in  partial  fractions.  The  result  is 


where 


- 1)  T 


2nABA  ffi  (Dj-C^)Oj^e  ^  erfci-a^^r) 


K  _  i) 

+  (D2  -  02)02  e  erfc  {-o^yfT ) 


+  erfc{~^) 
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(o^  -  1)  T 

-  O3  C3  e  eric  (-  0^^) 


(0=-l)T 

+  ^2  0^6  erfc  (o  ^  ^ 


Og  =  V  1-3  and  the  constants  are: 
1+0, 


C,  = 


C.  = 


c.  = 


c.  = 


D,  = 


Pr'2>K-‘?» 

O2-  (03-Oj)(03-l) 

1+  c:. 
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(A-15) 


1  (o^-02)(opl) 
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riiTii'iitiri  *•  ft*. 


^  — ;  wtvj  -  a,  . 


These  constants  have  the  following  useful  properties  which  can  be  verified 
by  direct  computation. 


S  C.  =  0 
i=l  ^ 


S  D.  =  0 
i=l  ^ 


Ij  a.  C.  =0 
i=l  ^  ^ 


£  0  D.  =0 
i=l  ^  ^ 


(A- 16) 


2  C.  =  1 
i=l  ^  ^ 


£  o^D.  =  1 
i=l  ^  ^ 


Also,  is  the  complex  conjugate  of  Cj  and  D2  is  the  complex  conjugate 
of  Dj^ . 

The  exact  solution  for  a  step  change  in  pressure  may  be  simply 
obtained  from  Equation  (A-14)  by  setting  C^,  C2,  and  all  equal  to 
zero.  Equation  (A-14)  is  not  in  the  form  best  suited  to  obtaining  numerical 
values  f<)r  m'/m,  for,  while  all  of  t'le  quantities  on  the  r.h.s.  of  the 
equation  are  complex  numbers,  m’/in  is  real. 

The  real  form  of  m'/m  may  be  found  as  follows.  The  definitions 
given  by  Equation  (4-62)  are  written  as 


Cj  =  x  +  iy 


=  x  -  iy 


(A- 17) 


and  the  following  identities  are  used  : 


er^c(-z)  =  1  -erfc(-z)  -  1  +  erf(z) 


A- 18) 


The  error  function,  itself,  can  be  expressed  in  series  form  as 


n  2n  + 1 
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-y/TT  ^_Q  n!(2n  +  1) 


(A-19) 


-f--  -r .  :>  ■"*- - ■-  '/ 


Using  Equations  (A-17  to  A-i9),  Equation  (A-14),  with  the  =  0,  can  then 
be  expressed,  after  much  algebra,  as 


iii:-2R  +  _i-ri+-~^  s  iriinn-ni 

2nABA  m  4AB  ^  n!{2n  +  l)  J 


Sr.  £ 


(A-20) 


n?.  ^j^^r<2"-^H[Rcos(2n+l)9-Ssl„(2n«)e] 


where; 


R  =  -  e^  (C  cos  vT  -  'n  sin  vt  ) 
S  =  -  e“  ^(ri  cos  vt  +  ^  sin  vt) 


u  =  -  1 


V  =  2  X  y 


ax-Sy  =  -§5g=a 
ri  =  Sx  +  ay 


(A-21) 


D  =  ,(x.-_U 
P  8yAB 


r  =  [-v-f+y^]^  =  (tc^o^)'^ 


9=  tan~^  y/x 


The  real  part  of  Equation  (A-14)  with  the  C.  0  (corresponding  to  the  solution 
for  an  exponential  pressure  change)  may  be  established  in  a  similar  manner. 
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The  value  of  the  response  function  obtained  from  Equation  (4-34) 
actually  corresponds  to  that  which  would  be  attained  after  many  cycles  of 
a  harmonic  disturbance.  During  the  period  immediately  following  the  initiation 
of  a  disturbance  the  value  of  the  response  function  may  be  quite  different 
as  a  result  of  transient  phenomena.  The  portion  of  the  response  resulting 
from  the  start  up  process  following  initiation  of  the  disturbance  damps  out 
exponentially,  and  the  limiting  value  given  by  Equation  (4-34)  is  asympto¬ 
tically  approached  as  time  increases.  The  number  of  wave  cycles  required 
for  these  initial  transients  to  damp  out  in  response  to  a  harmonic  disturbance 
varies  quite  significantly  as  a  function  of  A  and  B.  The  response  may  also 
initially  overshoot,  or  undershoot,  the  asymptotic  value. 

Examples  of  the  transient  burning  rate  response  to  a  pressure  disturbance 
of  the  form  AP  =  .  1  cos  nt  are  presented  in  Figures  A-1  to  A-3,  for  three  sets 
of  A  and  B  values: 

A  =  15,  B=  7 

A  =  11.5  B  =  .64 

A  =  20  B  =  .9 

The  »e  examples  illustrate  the  significant  effect  of  A  and  B  on  the 
initial  characteristics  of  the  response,  as  well  as  the  effect  on  the  ultimate 
magnitude  of  the  response  function. 
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